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1 Introduction 


In February of 1994, an effort from the Fluid Dynamics and Information Sciences 
Divisions at NASA Ames Research Center with McDonnel Douglas Aerospace 
Company and Stanford University was initiated to develop, demonstrate, vali- 
date and disseminate automated software for numerical aerodynamic simulation. 
The goal of the initiative was to develop a tri-discipline approach encompassing 
CFD, Intelligent Systems, and Automated Flow Feature Recognition to improve 
the utility of CFD in the design cycle. This approach would then be represented 
through an intelligent computational system which could accept an engineer’s 
definition of a problem and construct an optimal and reliable CFD solution. 

Stanford University’s role focused on developing technologies that advance 
visualization capabilities for analysis of CFD data, extract specific flow features 
useful for the design process, and compare CFD data with experimental data. 

During the years 1995-1997, Stanford University focused on developing tech- 
niques in the area of tensor visualization and flow feature extraction. Software li- 
braries were created enabling feature extraction and exploration of tensor fields. 
As a proof of concept, a prototype system called the Integrated Computational 
System (ICS) was developed to demonstrate CFD design cycle. 

The current research effort focuses on finding a quantitative comparison of 
general vector fields based on topological features. Since the method relies on 
topological information, grid matching and vector alignment is not needed in 
the comparison. This is often a problem with many data comparison techniques. 
In addition, since only topology based information is stored and compared for 
each field, there is a significant compression of information that enables large 
databases to be quickly searched. This report will (1) briefly review the tech- 
nologies developed during 1995-1997 (2) describe current technologies in the 
area of comparison techniques, (4) describe the theory of our new method re- 
searched during the grant year (5) summarize a few of the results and finally 
(6) discuss work within the last 6 months that are direct extensions from the 
grant. 

2 Review of Previous Work (1995-1997) 

This section reviews the various visualization techniques and tools developed 
under the grant during the years 1995 through 1997 that aid in the analysis of 
CFD data. Since this section is only a review, past papers are cited or enclosed 
that elaborate on the various techniques. 

The tools are general enough to be applied to any fluid data set and for the 
topology tools to a broader set of scientific data; however, for the purpose of 
this section one particular problem is focused upon: a hemispherical cylinder 
with an incidence angle of 19 degrees, an incoming Mach number of 1.2 and a 
Reynolds number of 445, 000. The extraction of shock waves and vortex features 
will be reviewed along with colored textures to simulate oil flow combined with 
pressure sensitive paint. In addition, a novel method for tensor visualization 
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using topological decomposition will be discussed. The section is completed 
by discussing how these tools are combined into an Integrated Computational 
System (ICS). 

2.1 Shock Wave Extraction 

A shock surface represents a sudden change of fluid properties. Typically, a 
shock is witnessed when a body travels at transonic speeds. The flow adjusts to 
the body by abruptly changing its direction with consequent changes in pressure, 
density, and temperature. This abruptness is caused by a nonlinear response of 
flow quantities such as pressure and density to the presence of a body or as a 
result of initial conditions. We have devised an effective method for shock ex- 
traction by applying knowledge of compressibility theory. By computing Mach 
numbers from velocities that are in the direction of the pressure gradient, and 
then extracting a corresponding isosurface, shock surfaces are identified. Noise 
and false detected surfaces (found passing through the boundary layer) are re- 
moved by enforcing the Rankine-Hugoniot relation, and by utilizing information 
about the geometry. A complete description of the algorithm is included in the 
supplemental entitled: Summary of Work on Shock Wave Feature Extraction 
in 3-D Datasets. The feature extractor was written using custom developed 
libraries that were specifically designed to ease the implementation of future 
extractors. This implementation reduced the solver’s code size, and provided 
management and visualization of the extracted feature. Figure 1 depicts the 
bow shock upstream of the hemispherical cylinder’s nose, along with a second 
shock on the top surface induced by the expansion due to the high angle of 
attack. 


2.2 Oil Flow Simulation 

3-D separated flows play a significant role in aerodynamics because of the close 
relationship between separation and the existence of vortices, and because of 
the effort placed in trying to avoid/understand them. Vortices are important 
structures of the flow far from the body as they affect lift and drag. The high- 
angle-of-attack flow about a hemisphere cylinder is a classical example of flows 
with massive separation. The first step in visualizing 3-D separated flows con- 
sists of depicting the structure of the vector field near the body, for which an 
adequate description is inferred from the skin-friction field — i.e., the 2-D tan- 
gential velocity field one grid plane away from the surface of the body. This 
method mimics a visualization method that is commonly used by experimental- 
ists, namely, the oil flow method. 

A recent major improvement was made by adopting a new rendering tech- 
nique, namely anisotropic textures. By generating textures, we render 2-D 
streamlines directly, without the need for integrating the vector field (Refer- 
ence [1]). Textures provide an adequate solution to problems of accuracy, do 
not require an integration of the whole vector field, and create continuous im- 
ages. 
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Figure 1: Shock extraction about a hemisphere cylinder. 


By combining coloring based on pressure values on the body surface this 
technique of oil flow simulation provides the researcher with an image that 
can be analyzed using critical point theory and can be compared with images 
from experiments. Figure 2 contains such a description of the flow. From the 
image one can clearly identify the lines of separation and re-attachment along 
with the critical points on them. Details of the algorithm are described in the 
supplemental entitled: Flow Visualization with Textures. 

2.3 Vortex Extraction 

Vortex identification is an important flow feature in the design of aircrafts. 
Locating vortices can aid in improving aircraft safety, and airport landing and 
take-off throughput. It can also help increase turbo-machinery efficiencies and 
reduce airframe vibration and audible noise. 

One method for vortex detection is to use topological skeletons. Separation 
fines can be located using critical point theory. These surfaces of separation and 
the associated vortices are extensions in the flow of the skin-friction topology. 
The surfaces of separation emanate from the body along lines of separation. 
They are computed by advecting a front of streamlines initially positioned along 
the lines of separation. Figure 3 shows the surfaces of separation which roll up 
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Figure 2: Oil flow simulation combined with pressure sensitive paint 


to form the vortices. 

2.3.1 Vortex Core Extraction 

There are several prevalent techniques for the identification of vortex cores. 
Critical point theory [2], pressure/density minima and helicity are just a few 
techniques used to extract vortex cores. Observations relating the medium 
eigenvector of the deformation tensor to the vorticity vector [3, 4] have lead 
to investigations into techniques based on the topology of tensor fields and our 
idea of hyperstreamlines. Our method for vortex core extraction is based on 
using the helicity density gradient [5] . This method is effective in subsonic flows 
where density methods typically fail. Unlike methods that use enstrophy, the 
sign is preserved hence allowing the identification of secondary vortices via the 
swirl direction. This method also produces a continuous feature by integrating 
along the helicity density gradient vector, unlike methods based on critical point 
theory which use interpolation. Figure 4 depicts one of the four vortex cores 
(in blue), and the relation to the surrounding flow via hyperstreamlines of the 
reversible momentum tensor. 
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Figure 3: Stream surfaces depicting separation topology 


2.4 Tensor Visualization 

Tensor data sets are at the heart of many engineering and physics disciplines, 
yet few methods have been devised for understanding and visualizing such data. 
In particular, second-order tensor fields are central to fluid mechanics. For 
example stresses, strain rate, and Reynolds stresses are all tensor quantities. 
In traditional approaches, selected data are displayed using simple local icons 
depicting, for example, eigenvalues and eigenvectors at selected points. Through 
this data reduction process, valuable information is lost or discarded. 

Since a tensor field is a continuous field, it should be visualized as such. The 
approach we have taken is to trace the trajectories of the eigenvectors of the 
tensor field (referred to as hyperstreamlines [6]). Any symmetric tensor data set 
can be represented locally by a set of orthogonal eigenvectors and eigenvalues 
which contain all the information about the tensor data at a given point in space. 
Hyperstreamlines are tangent curves to the principal eigenvectors; the principal 
eigenvectors are ordered according to the magnitude of their eigenvalues. After 
a particular eigenvector is chosen, integrating along the direction of the chosen 
eigenvector, a tangent curve is generated along which we can display, in an 
orthogonal manner, the remaining two eigenvectors. This can be done as a set 
of orthogonal vectors or ellipses having principal axes in proportion to the two 
eigenvalues. 
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Figure 4: Vortex core (solid blue) encompassed by hyperstreamlines (major 
eigenvector) of the reversible momentum tensor 



In this manner, a display is generated that takes advantage of the continu- 
ity of the data, and all information in the original data are retained. Figure 5 
depicts the reversible momentum tensor in a clear and concise form. Hot col- 
ored regions depict kinetic energy density, while the hyperstreamlines trace the 
flow’s path. Details of using tensor visualization methods in the CFD process 
are discussed in the following enclosed supplements: The Topology of Symmet- 
ric Tensor Fields [7], The Topology of Three-Dimensional Symmetric Tensor 
Fields [8], The Topology of Symmetric, Second-Order 3-D Tensor Fields [9]. 

2.5 Integrated Computational System (ICS) 

Over the past few years, the development of the techniques described above 
were largely funded by NASA Ames with additional NSF funding for the work 
in visualizing the topology of 3-D tensor fields. To unify the visualization effort 
based on the principle “Analyze then Visualize”, a framework for solving flow 
features was designed recently using an object-oriented methodology. Exploiting 
encapsulation, inheritance, and dynamic function binding, we have built a com- 
prehensive and extendible system, known as Integrated Computational System 
(ICS). The tools in ICS have been applied to solving problems in computational 
fluid and solid mechanics and can be used to compare results of multiple flow 
solutions. Schematically this system is shown in Figure 6. 

Each feature extraction tool is a separate inherited solver. Each feature 
solver inherits data and methods from a base feature solver, and are all managed 
by a feature solver manager. Problems are submitted to the manager in a 
generalized, typically attribute oriented way. The feature solver manager makes 
the decision and returns a specific feature solver that is best suited for the 
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Figure 5: Rake of hyperstreamlines of the major eigenvector field of the re- 
versible momentum tensor. 


problem description. This returned solver being a child of a base feature solver 
has a common interface. This provides a great benefit in terms of minimizing 
documentation, and standardizing the interaction with other feature solvers that 
will be created. This polymorphic behavior has been extended to solutions as 
well. Each feature solver returns a solution which is a feature that has inherited 
from a base feature class, known as a primitive. A standard interface exists for a 
primitive to display itself to the screen, or in other formats such as VRML 1 . The 
primitive is typically small, and usually contains some compression such that 
the transfer to the client’s machine requires a minimal amount of bandwidth. It 
is important to note that a feature solver manager provides a division between 
the interface and the implementation. The manager can be running on the 
same system as the client code or on another computational server in a different 
connected network. 

The prototype currently provides the means to conduct a simulation using a 
full Navier Stokes flow solver, and allows automatic extraction of flow features 
such as shocks, vortex cores, and simulated oil flow combined with pressure sen- 
sitive paint. It also includes a tensor visualization module that creates tensor 
quantities relevant to fluid mechanics such as the velocity gradient, deformation, 

Virtual Reality Modeling Language provides a standard for three-dimensional data to be 
shared across the internet 
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Integrated Computational System (ICS) 


3-D Navier Stokes Solver 
Visualization tools 
Graphical User Interface 
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>Boussinesq problem 
oStresses in weld 


Figure 6: Visualization tools developed by our group at Stanford University. 


and viscous stress tensors. In addition, user defined tensors can be loaded and 
explored within the system. All data is visualized using a custom developed 
application known as FASTlook which is similar to Silicon Graphic’s Open In- 
ventor’s ivview. In contrast, however, FASTlook is designed for visualization of 
data-rich primitives found in scientific visualization. In addition, at the request 
of NASA Ames, FASTlook has been converted to a network helper application 
with the intention of incorporating the software within web based applications 
such as DARWIN [10]. 

Supplemental documentation is included that describes FASTlook entitled 
Pnm View: A Scientific Visualizer User’s Guide c2.0 and documentation de- 
scribing the use of the solution primitives entitled PrimView: A Scientific Vi- 
sualizer Library Guide v2.0. 


3 Comparison Techniques 

There exist a variety of comparison techniques for vector fields. These tech- 
niques basically fall into three general categories: Image, data, and feature 
extraction based comparisons. In most of these cases, comparisons are made 
visually [11]. Image based comparisons work on the computer generated image. 
Often times, a numerical data set is converted into an image that simulates an 
experimental visualization technique (computational flow imaging). This may 
be easier than extracting a vector field from an image, such as Schlieren. How- 
ever, visualizing a field in 3-D is quite difficult. Often times, these techniques 
are limited to two dimensions. In addition to side-by-side comparison of images, 
other techniques include image fusion, and Fourier analysis [12]. 

Data level comparison techniques operate directly on the raw data. An accu- 
rate comparison requires proper grid alignment which can involve problematic 
interpolation between two fields [13]. 

The last comparison category is the extraction of features . Typically fea- 
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tures axe flow specific such as vortex cores, shock surfaces, or topology. Often 
times there is a geometric representation of the feature and possibly a semantic 
representation of the system which can be compared using a pattern recognition 
technique [14]. This may lead to more robust comparisons. Past study in our 
group has focused on the geometric structure of vector fields [15]. However, this 
geometric structure can be visually deceiving since two vector fields may have 
the same underlying topological structure but are dissimilar in appearance [16]. 
Therefore, a quantitative measurement for comparison of vector fields is essen- 
tial. 

4 Description of a Vector Field 

A 2-D vector field can be described as a system of two simultaneous differential 
equations having the following form: 


where F and G are continuous and have continuous partial derivatives in some 
region D. 

A vector field is typically described by the number, type, and arrangement 
of critical points (or equilibrium points). These points are where the system is 
defined to be F(x, y ) = 0, G(x, y) = 0. The number and nature of critical points 
will not change under continuous transformation. A critical point is said to be 
isolated or simple if there is an open neighborhood around it that contains no 
other critical points. For this report, we focus entirely on simple critical points. 
The global topology of the vector field is defined as the critical points and the 
set of their connecting streamlines. These streamlines (separatrices) divide the 
field into regions that are topologically equivalent to uniform flow. Hence, only 
the topology is needed to reconstruct the field and therefore is useful as a means 
of differentiating vector fields. 

4.1 Classification of Simple Critical Points 

The behavior of the flow about a critical point can be analyzed by investigating 
the streamlines in the neighborhood of the critical point. If we are sufficiently 
close to the critical point (say a distance dx,dy away) in most cases a first order 
Taylor series expansion of the velocity field is sufficient: 


Vx = ~dt = F ( x ’ y ' ) 


( 1 ) 


v v = ^ = G ( x ^y) 



(2) 
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Attracting Ftcus 
a < 0 
0 < 0 


Center 
a =0 
0 < 0 




Attracting Star 
a <0 
p =0 


A t t ra cti n g Neds 
a < 0 
p >0 


Figure 7: basic patterns for simple critical points 


Hence, the flow pattern is completely determined by the 2X2 Jacobian ma- 
trix, Jij — |2i. = i 2) evaluated at the critical point location. The various 

patterns formed in the phase-plane space can be seen by analyzing the eigenval- 
ues of the Jacobian. The patterns are sketched in Figure 7. Notice a positive or 
negative real part (denoted by a) is indicative of repelling/ attracting behavior. 
And if an eigenvalue has an imaginary part (/3 < 0), it indicates circulation 
about the point, otherwise asymptotic behavior is exhibited. 


5 Vector Field Representation 
using Clifford Algebra 

In [17] [18], Sheuermann et al. introduced Clifford algebra for vector field visu- 
alization. Clifford algebra provides a nice way to describe the relation between 
real and complex numbers in 2D space. The vector fields are defined over a 
complex field in this algebra and the nonlinear vector fields are represented as 
multiplications of linear fields. 

For the Euclidean plane we get a 4-dimensional R-algebra (?2 with the basis 
ljei>e 2 ,i = eie 2 as a real vector space. Multiplication is defined as associative, 
bilinear and by the equations 


le, 

CN 

T— 1 

II 

II 

(3) 

eiej 

= U = l,2 

(4) 

eie 2 

- e 2 ex 

(5) 

l 2 

i-H 

II 

H 

II 

(6) 

A 

= W = 1,2 

(7) 
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-1 


(8) 


i 2 = 


with 



The usual vectors ( x , y) £ R 2 are identified with 

xe\ + ye 2 £ E 2 C G 2 (9) 

and the complex numbers a + bi £ C with 


al + bi £ G 2 


(10) 


5.1 Vector fields in Clifford space 

A Clifford vector field is just a multivector field with values in R 2 C Go 

v : R 2 -> R 2 C G 2 (11) 

Let z = x + iy, z = x — iy be complex numbers in the Clifford algebra. This 
means 

x = -( z + z) (12) 

y = ^r(z-«) ( 13 ) 

We get 

v(r) = v 1 (x,y)e 1 +v 2 (x,y)e 2 

(l( z + z),±(z-z)^ ei 

= E(z,z)e 1 (14) 

Generally, a linear vector field can easily be shown as: 
v(r) = E(z, z)e 1 

= ( az + bz + c)e 1 (15) 


where a,b,cC C. 

Let E : C 2 — * C be the polynomial so that v = E(z,z)e\. Let F* : C 2 -> 
C,k = 1, . . . , n be the irreducible components of E, so that E(z, z) — rifc=i 
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Table 1 : Classification of Critical Points using a B values 


a 

B 

Type 

a 

B 

Type 

= 0 

< 0 

Center 

— 

> 0 

Saddle \ B \ > |a| 

> 0 

< 0 

Repelling Focus 

< 0 

< 0 

Attracting Focus 

> 0 

= 0 

Repelling Star 

< 0 

= 0 

Attracting Star 

> 0 

> 0 

Repelling Node \ B \ < |a| 

< 0 

> 0 

Attracting Node \ B \ < |a| 


then an arbitrary polynomial vector field with isolated critical points can be 
expressed as: 


v(t) = E{z,z)e i 

n 

= JJ(a*z + f>*z + c*)ei ( 16 ) 

*=i 

where z* is the unique zero of a^z + b^z 4 - c* . 


6 a-fl Space and its Use as a Metric 

For a linear vector field v = (az + bz + c)e i, let a = ai 4 - a-ii and b ~ b\ + b^i. 
Eigenvalues of the Jacobian around its critical point zo are Aj — b\ + y/\a\ 2 — b\ 
and A2 = 61 — \/M 2 — b$- 

Let a = 61 and B = sign(\a \ 2 — 6|)v / ll a l 2 — &il> criteria for basic patterns of 
simple critical points axe: 

Selection of a and B as shown in Figure 7 and delineated below can be 
mapped to ai, 61, 02, 62 to yield any desired field: 

Notice our definition of saddle is more relaxed than shown in Figure 7 . The 
values of a and B determine the type of critical point but it is not sufficient to 
be used as a metric to differentiate between two types of critical points. So we 
introduce a new a~B space where the 8 simple critical points are mapped onto 
the a, B axes at their respective (a, B) points. Vectors in this space obey all 
the rules defined for a regular 2 -D Euclidean space. All points in this space are 
normalized as follows: 


a = 


y/a 2 +0 2 


B' = 


y/a 2 +0 2 


( 17 ) 


It is shown in [ 19 ] that the actual values of a and B do not determine 
the portrait of the critical point only the ratio between them. Hence, this 
normalization maps all points onto a unit circle (Figure 8) and thereby provides 
a means of relatively quantifying the difference between various points. Also 
note that a regular vector field with no critical points, v = const ■ e\ has a = 0 
and B = 0 and sits at the origin of the unit circle. For the remainder of the 
report, a and B values will be assumed normalized. 

A multiple point with a set of a’s and / 3 ’s corresponds to a set of points in 
the a — B space. For example, v = z 2 e 1 is a dipole which has two (1, 0) point 
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Figure 8: Basic patterns for critical points in a- 8 space; C for center, RN for 
node, AN for attracting node, RF for repelling focus, AF for attracting focus, 
St for star, Sa for Saddle and R for regular point 


in a — 8 space; and v = (z — (2 + 2 i)z + ci)(z + (2 + 2 i)z + C 2 )ei has one point 
at ( — ^|) and another point at (-^-,-^|) hi a — 8 space. 

7 Earth Mover’s Distance 


7.1 EMD analysis 

The Earth Mover’s Distance is first introduced in [20] [21] for content-based 
image retrieval in a large data base. It is used to compute the minimal amount 
of work that must be performed to transform one feature distribution into the 
other. Feature distribution in [20] [21] are the color and texture signatures of 
an image. 

After careful study, we found that the EMD concept can be used to compute 
the differences between vector fields. Here, the feature distribution is redefined 
as the characteristics of a vector field. 


Definition 1 (feature distribution) A feature distribution for a vector field 
is the set of a and 8 values associated with the vector field’s critical points: 

{(Qi , 8i), (Q2, 82), ■ ■ ■ , (an, /?„)}. 

Definition 2 (Energy) The energy for a vector field is: 


Energy = 


n 


N 


!>?+©• 
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where n is the total number of critical points in this field. 

This energy here is a quantity that characterizes the critical points of a vector 
field. It is different from the physical energy. The concept ’’work” is used to 
measure the energy differences between two vector fields or the amount of energy 
used to transform one vector field into the other. 


Definition 3 (Work) For two vector fields with feature distributions 


and 




The amount of work nec essary for transforming one ve ctor field into the other 
is defined as: Work = VH!!=i(( a > - a') 2 + (0i~ P[) 2 )- 


Intuitively, given two feature distributions, one distribution can be seen as a 
set of discrete point-objects with a certain amount of mass of earth spread in 
space, the other as a collection of holes in the same space. The work measures 
the least amount of energy needed to fill the holes with earth and is called the 
Earth Mover’s Distance (EMD). Computing the EMD is based on a solution to 
the old transportation problem from linear optimization [22]. This is a bipartite 
network flow problem which can be formalized as the following linear program- 
ming problem: Let I be a set of suppliers, J a set of consumers, cy the cost to 
ship a unit of supply from i £ I to j £ J 


Cij = yj (a i - clj ) 1 + (ft - 


and it is the same as the Euclidean distance dij = ||iT; — Vj\\ in a- (3 space. A 
critical point either exists as a whole or does not exist, it can not be split. In 
this case, the transportation problem has the property that the optimal flow /y 
can only be 0 or 1 [23] . We want to seek a set of /y that minimizes the overall 
cost: 

EMDfx, y) = miri^T ^ Cy/y (18) 

iel jeJ 

subject to the following constraints: 


fa > o iel, j £ J 

(19) 

^ ^ fij = Vi i j G J 

(20) 

i€l 


= Xi, i£l 

(21) 

jeJ 


= J2 Xi 

(22) 

jeJ iel 
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Where Xi is the total supply of supplier i and yj is the total capacity of consumer 
j. Constraint (19) allows shipping of supplies from a supplier to a consumer 
and not vice versa. Constraint (20) forces the consumers to fill up all of their 
capacities and constraint (21) limits the supply that a supplier can send as a 
total amount. Constraint ( 22) is a feasibility condition that ensures that the 
total demand equals the total supply, in other words, the distributions have the 
same overall mass and the EMD is a true metric [20]. 

It is likely that a set of vector fields will not have the same number of 
distributions. In order to satisfy constraint (22), we can create regular points 
to make the supply equal the demand without changing the vector fields. For 
example, if the supplier field contains 3 critical points 

3 

v = JJ(aiZ + + Cj)ei (23) 

i=l 

and the consumer field contains 5 critical points 

5 

v ' = II + b 'i* + C i) e l ( 24 ) 

j=l 

The supplier side has two fewer points in the a — 0 space. Now let 

3 

v = JJ(a,z + biz 4- a) ■ 1 • lei (25) 

i = 1 

and the vector field remains unchanged. However, now we have two more regular 
points corresponding to 1 with a = 0 and 3 = 0, and both the supplier and 
the consumer have 5 points in their feature distributions. All the conditions are 
satisfied, and we are ready to compute the EMD for these two fields and find 
out the dissimilarity between them. 

In order to evaluate the meaningfulness of our new metric, we use Multidi- 
mensional 5caZin^(MDS) [24] [25] to embed the vector fields in a two-dimensional 
Euclidean space so that distances in the embedding are as close as possible to the 
true EMDs between vector fields. The MDS is introduced in the next section. 

8 Display of EMDs for a Large Set of 
Vector Fields 

The above discussions are for comparison of a pair of vector fields. If there 
exist a large set of vector fields and we want to compare their topologies, it is 
necessary to display them in a more meaningful way than a ID list sorted by 
their EMDs. Yossi Rubner et al. have used Multidimensional Scaling Method 
(MDS) [20, 21] to display a set of images on a 2D map. Given n objects in a 
high dimension, the MDS method computes a configuration in a lower dimension 
space such that the distance between every pair of objects in this low dimension 
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Figure 9: EMD of the flow over an a) airfoil b)circular cylinder. 


space best matches the real distance in the high dimension. Inspired by their 
work, we compute the EMDs between every pair of vector fields and position 
the vector fields on a 2D map such that the distances between the vector fields 
match their EMD values as close as possible. 


9 Application: Flow over an Airfoil and Cylin- 
der 

Rogers and Kwak computed the flow past a 2-D airfoil at —90° angle of at- 
tack [26]. The model was of interest since the flow of the wake of an XV-15 Tilt 
Rotor aircraft degraded the lifting capability during hover. An incompressible, 
time accurate, Navier-Stokes code with artificial compressibility at a Reynolds 
number of 200 was used to compute the flow over a NACA 64A223M airfoil. 
Fifty frames were computed. During this time the flow entered into a peri- 
odic vortex shedding cycle. Earth mover’s distance was computed over the 50 
frames. The plot in Figure 9a depicts the EMD comparison of frame 1 with 
the remaining 49 frames. At frame 1, the EMD is zero, which is expected since 
the work required to convert a frame into itself is zero. The periodic nature 
is apparent. We see a repetition approximately every 17 frames. Also we see 
a sudden EMD rise when comparing frame 1 with frame 8 indicating a signif- 
icant topological feature difference. Frame 1 contains three critical points: an 
attracting/repelling focus and a saddle. Frame 8 contains 5 critical points: two 
saddles, an attracting/repelling focus pair and a node. Since the flow is incom- 
pressible, the velocity divergence, V • v, is expected to be zero everywhere in 
the flow. Hence only saddles and centers are to be extracted. However, it is 
common due to numerical computation for a to not exactly be zero, however, 
we should not expect to find a node. Upon closer examination of the data, the 
velocity divergence for certain frames is not zero near the tips of the foil where 
nodes are being extracted. We believe that the flow solver may not have fully 
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Figure 10: Topologically similar frames 1, 34 and 50 of flow about an airfoil. 
Frame 8 is topologically dissimilar. 


converged and therefore we see this sudden jump of discontinuity in the field. 
The LIC images of frames 1, 34, and 50 are depicted in Figure 10 have very 
similar earth mover’s distance and as can be seen look nearly identical. Frame 
8 differs from the others due to its variation in topology (formation of nodes) 
and is apparent in the figure. 

We contrast the flow over an airfoil with the flow over a circular cylinder 
simulated by Rogers and Kwak under the same flow conditions [26]. In this case, 
the flow is divergence free and the EMD values are quite similar. Thirty frames 
were computed capturing a complete cycle of vortex shedding. As can be seen 
from the plot in figure 9b, Frames 1 and 16 have nearly identical EMD values 
leading one to believe the period to be every 15 frames. Due to the symmetry of 
the flow, this is not far from the truth. In fact the flow produces a mirror image 
of itself every 15 frames as it sheds the alternate vortex and hence leads to the 
same topology. Figure lib depicts the alternate vortex being shed to the image 
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Figure 11: Topologically similar frames 1 and 16. FYame 4 depicts the down 
stream dissipation of saddle-center pair producing a larger EMD. Frame 6 de- 
picts formation of new saddle-center pair at the bottom of the cylinder. 


found in figure 11a. Furthermore, we see from figure 9b an increase in EMD 
value for frame 3. This increase is due to the dissipation of the saddle-center 
pair as it moves down stream (figure 11c). The EMD drops by frame 6 as the 
next saddle-center pair is shed (figure lid). By frame 16, the saddle-center has 
moved down stream such that the a, (3 values are nearly identical to frame 1. 
Two frames later the saddle-center dissipate and the cycle repeats. 


10 Discussion 

We have demonstrated the effectiveness of topology based feature comparisons 
for vector fields. The use of a quantitative measure between fields provides the 
means for fast automated comparisons as well as an indepth study of flow fields 
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as demonstrated with the time history data. For the airfoil data, we have shown 
the effectiveness of the method as a diagnostic tool. The clear EMD difference 
provides an immediate alert into calculation problems for particular frames. For 
the cylindrical data, the periodic nature of the flow was revealed. The EMD 
difference also provides insight into the evolution of the flow field. 

11 Direct Extensions from the Grant Research 

The work conducted under the grant introduced a quantitative method for com- 
paring 2-D vector fields based on the number and type of critical points that 
comprised the field. However, the arrangement of the critical points was not 
considered, potentially causing two very different fields with the same type of 
critical points but different streamline connections to be detected as similar 
fields. The work has been extended by considering the connections between 
critical points, thereby improving the representation of the vector field. The 
vector field’s topology is represented using an attributed relational graph, and 
with the use of conventional graph matching algorithms a comparison of the 
fields can be made. The supplemental paper entitled Topology Based Vector 
Field Comparisons Using Graph Methods [27] discusses these graph methods. 

Furthermore, the work under the grant focused on 2-D simple critical points. 
Many problems occur in three dimensions, and therefore a technique for 3-D 
comparisons has been recently developed. The extension to three-dimensions 
follows the path of our previous work, rethinking the representation of a critical 
point signature and the distance measure between the points. The discussion 
of this technique is included in the supplemental entitled Feature Comparisons 
Of 3-D Vector Fields Using Earth Mover’s Distance [28]. 
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Abstract 


Scientific visualization transforms complex data obtained from experiments, obser- 
vations and numerical simulations into an image amenable to understanding by the 
human visual system while maintaining the integrity of the information. It is a means 
to communicate between human perception and abstract physical world. 

Vector (first-order tensor) and second-order tensor fields are multivariate, mul- 
tidimensional data sets that are typically very difficult to comprehend. Visualizing 
such complex data is a very important subfield in scientific visualization. 

Much research has been accomplished on vector field visualization. However, 
virtually no work has been carried out on quantitative comparisons of similarities 
and differences between vector fields. This dissertation introduces a novel approach 
to define a topology based measurement for such a purpose. The usefulness of this 
measurement can be seen when comparing computational and experimental flow fields 
under the same conditions. Furthermore, its applicability can be extended to such 
cumbersome tasks as navigating through a large data base by searching for a similar 
topology. 

This new.measure relies .on .the. .use of. critical points, which are a key feature of 
vector field topology. In order to characterize critical points, a and (3 parameters are 
introduced. They are used to form a closed set of eight unique patterns for simple 
critical points. These patterns are also basic building blocks for higher order nonlinear 
vector fields. In order to study and compare a given set of vector fields, a measure of 
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distance between different patterns of critical points is introduced. The basic patterns 
of critical points are mapped onto a unit circle in a - (3 space. The concept of Earth 
Mover’s Distance is used to compute the closeness between various pairs of vector 
fields, and a nearest-neighbor query is thus produced to illustrate the relationship 
between the given set of vector fields. 

Very few methods have been developed to understand and visualize second-order 
tensor fields due to their complex nature and large data sizes. Topology based study 
has been done in two-dimensional space. This dissertation extends the research into 
three dimensions, an area that hasn’t been studied before to the best of our knowledge. 

The basic constituents of second-order tensor fields are degenerate points. They 
play a role similar to critical points in vector topology. In this dissertation, we address 
the conditions for the existence of degenerate points and based on these conditions 
we predict the distribution of degenerate points inside the field. Every tensor can 
be decomposed into a deviator and an isotropic tensor. A deviator determines the 
properties of a tensor field, while the isotropic part provides a uniform bias. De- 
viators can be three-dimensional or locally two-dimensional. The triple degenerate 
points of a tensor field are associated with the singular points of its deviator and the 
double degenerate points of a tensor field have singular local 2-D deviators. Con- 
trol functions are in charge of the occurrences of a singularity of a deviator. These 
singularities can further be linked to important physical properties of the underlying 
physical phenomena. For example, for a deformation tensor in a stationary flow, the 
singularities of its deviator actually represent the area of the vortex core in the field; 
for a stress tensor, the singularities represent the area with no stress; for a viscous 
flow, removing the large, isotropic pressure contribution enhances dramatically the 
anisotropy due to viscosity. 
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Chapter 1 


Introduction 

1.1 Scientific Visualization 

“The purpose of computing is insight, not numbers” [1] 

[Richard Hamming, 1962] 

Scientific visualization transforms complex data obtained from numerical simula- 
tions, observations and experiments into a graphical abstraction and renders it into 
an image amenable to understanding by the human visual system while maintaining 
the integrity of the information. 

The concept of visualizing numerical data did not appear recently. From the 
middle of 17th century' to the beginning of 20th century, there have been techniques 
developed by some of the greatest scientists including Halley, Watt, Descartes, Lam- 
bert,. Playfair, and. von Humboldt. -However, it was -not until the introduction of 
computer graphics in the 1960s, when data visualization was brought into the fore- 
front of science and engineering. In particular, McCormick’s report on Visualization 
in Scientific Computing [2] in 1987 stimulated an explosion of research in this “new” 
field. Today, computer generated images have become a standard for communication 
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involving human perception and abstract representation of the physical world. 

Multidimensional Multivariate Visualization 

Multidimensional multivariate (mdmv) visualization is an important subfield of sci- 
entific visualization. It was studied separately by statisticians and psychologists long 
before computer science was deemed a discipline [3]. This subfield studies multi- 
ple parameters and the key relationship between them. Multidimensional refers to 
the dimensionality of the independent variables, and multivariate refers to the di- 
mensionality of the dependent variables [4]. For example, a deformation tensor field 
(symmetric tensor) observed and recorded in a three-dimensional space at various 
locations produces 3d6v data. 

Scientists have studied multivariate visualization since 1782 when Crome used 
point symbols to show the geographical distribution in Europe of 56 commodities [5]. 
However, before Tukey’s exploratory data analysis, tools for multivariate visualization 
usually consisted of colored pencils and graph paper, and mainly dealt with small- 
sized one or two variate data. The appearance of low-priced personal computers and 
workstations during the 1980s breathed new life into graphical analysis of mdmv data. 
During the last decade, hundreds of new mdmv visualization techniques have been 
invented [6]. Among them, a large number have been devoted to vector and tensor 
field visualization such as described in [7, 8, 9, 10, 11, 12, 13, 14, 15]. 

1.2 Vector and Second-order Tensor Fields 

Vectors and second-order tensors are both multidimensional, multivariate quantities. 
A vector field in a three-dimensional space has three components at each location 
in the field, therefore generates 3 d3v data; whereas a general tensor field in a three- 
dimensional space has nine components and generates 3d9v data. In some special 
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cases, the number of independent components of a tensor is reduced. For example, a 
symmetric tensor has only six independent components and an anti-symmetric tensor 
has only three, therefore produce 3 d6v and 3d3v data, respectively. 


1.2.1 Vectors, Second-order Tensors and the Field Concept 

A tensor is a quantity whose key property is the transformation law of its components, 
i.e, the way its components in one coordinate system are related to its components in 
another. The order of a tensor describes its complexity [16]. The simplest tensors are 
a scalar and a vector which are a zeroth-order and a first-order tensor, respectively. 
A second-order tensor is next in order of complexity in the tensor family. The focus 
of this dissertation will be on vectors and second-order tensors. 1 


Vector 

A vector is a quantity uniquely specified in any coordinate system by three real 
numbers (or three components). The components of a vector transform under changes 
of the coordinate system according to the law (Einstein summation rule is used here): 

V! = a ik V k 

where Vf, V k are the components of the vector in the old and new coordinate systems 
K and K’, respectively, and is the cosine. of.the .angle .between the ith axis of K’ 
and the kth axis of K. Generally, if a vector vanishes in one coordinate system it 
vanishes in other coordinate systems as well [16]. 

: in this dissertation, “second-order tensor” will be very often referred as just “tensor” for 
convenience. 
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Second-order Tensor 

A second-order tensor is a quantity uniquely specified by nine real numbers (the 
components of the tensor) which transform under changes of the coordinate system 
according to the law: 


T'ik — CXil&kmTim 


( 1 . 1 ) 


where Ti m and are the components of the tensor in the old and new coordinate 
systems K and K’, respectively, and an is the cosine of the angle between the ith axis 
of K’ and 1th axis of K (similarly for a*m) [16]. Generally, if all the components of a 
tensor vanish in one coordinate, they also vanish in any other coordinate systems. 

A second-order tensor is often expressed as a matrix: 


^ r P *~P r T' ^ 

-'ll -t 12 I 13 

T 2 \ t 22 t 23 

^ ^31 T 32 T 33 y 


( 1 . 2 ) 


Tensor Field Concept 

A tensor field is a rule assigning a unique value of a tensor to each point of a certain 
volume V in space. Let r be the position vector of a variable point of V with respect 
to the origin of some coordinate system. Then a tensor field of order n is indicated 
by 

= ( r ) ( 1 - 3 ) 


Thus, a vector field is defined as: 


T it = V h (r) 


(1.4) 
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And a second-order tensor field is defined as: 


T 


W2 


( T n (x, y, z) T 12 (x, y, z) T 13 (x, y, z) ' 

T 21 (x, y, z) T 22 (x,y,z) T 23 (x, y, z) 

K T 3 i (x, y, z) T 32 (x, y, z) T 33 (x, y, z) ) 


(1.5) 


Definition 1 (Critical Points) A critical point is a point in a vector field where 
all three vector components are zero and the streamline slope is indeterminate. 


Critical points are the only points in a vector field where tangent curves may 
cross each other. In the case of a steady state velocity field, tangent curves represent 
streamlines. Critical points are characterized according to the behavior of nearby 
tangent curves. There are eight basic patterns for simple (first order) critical points 1 , 
namely, an attracting/repelling star, an attracting/repelling node, a saddle, an at- 
tracting/repelling focus and a center. A particular set of tangent curves-separatrices 1 , 
which end on critical points, are of special interest because they define the skeleton 
which characterizes the global behavior of all other tangent curves in the vector field. 
The number and nature of critical points in a vector field remain unchanged under a 
continuous transformation. 


Definition 2 (Degenerate Points) A degenerate point of a tensor field T : E — > 
C ( R m , R m ), where E is an open subset of R m , is a point xo € E where at least two 
of the m eigenvalues of T are equal to each other [38]. 

In the case of two-dimensional tensor fields, there are only two eigenvalues A x 
and A 2 , and x 0 is a degenerate point iff Ai (x 0 ) = A 2 (x 0 ). For three-dimensional 
tensor fields, various types of degenerate points exist, corresponding to the conditions 
Ai (x 0 ) = ^2 (xo), ^2 (x 0 ) = A 3 (x 0 ), or Ai (x 0 ) = A 2 (x 0 ) = A 3 (x 0 ). 


1 will be explained in Chapter 2 
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1.2.2 Combinatorial Topology of Vector and Tensor Fields 

Combinatorial topology', also known as “rubber sheet geometry”, is a branch of ge- 
ometry. It studies the properties of figures that endure when the figures are subjected 
to continuous transformations. A topological property of a figure is a property pos- 
sessed alike by the figure and all its topological equivalents. In combinatorial topology, 
complicated figures are constructed from simple ones and their properties are deduced 
from the simple figures [17]. 

Combinatorial topology has extensive applications in geometry and analysis, many 
of which result from connections with the theory of differential equations. It may 
seem surprising that such superficially different subjects as topology and differential 
equations could be related, but research has shown that a link between these two 
subjects is the concept of a vector field [17]. Recent developments in scientific visual- 
ization have shown that vector fields and their topological structures also play a very 
important role in analyzing second-order tensor fields. 

Second-order tensors are fully represented by their eigenvectors and associated 
eigenvalues. 

Tel = A id (1.6) 

where A, and e) (z = 1,2,3) are eigenvalues and eigenvectors of the tensor T, respec- 
tively. The Aj’s represent all the amplitude information while the efs represent all 
the directional information of T. Visualizing a tensor field is equivalent to visualizing 
its eigenvector fields. However, unlike a vector field, eigenvectors are vectors with 
sign indeterminacy._This.remarkableJeature.distinguishes eigenvector fields from or- 
dinary vector fields and makes their topological features even simpler. The basic 
constituents of tensor topology are degenerate points where at least two of the eigen- 
values are equal to each other. They play a role similar to critical points in vector 
fields. 
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1.3 Motivations and Objectives 

1.3.1 Overview of vector field visualization 

Vector fields have vast applications in physical sciences. The force fields arising from 
gravitation and electromagnetism; the velocity vectors of a fluid motion, found for 
example over an airfoil; and gradients, such as the pressure gradient appearing on a 
weather map are all elements of vector fields. 

Vector field visualization has received an enormous amount of attention. Tra- 
ditionally, visualization has been part of measurements in complex physical experi- 
ments. Experimental fluid mechanics relies heavily on physical visualization for flow 
measurements and analysis. Many techniques have been developed in the past: seed- 
ing a flow with smoke or dye dates back decades [18]; using hundreds of tiny tufts 
attached to screens for visualization of two-dimensional cross sections of flows [19]; 
and photographs of oil streak patterns on body surfaces. These techniques provide 
a wealth of information for both local and global flow structures and have had a 
profound influence on flow modeling studies. However, direct results from experi- 
mental visualization usually give us only a qualitative measure and the derivation of 
quantitative information is accomplished by image processing and analysis. And more 
recently, with the fast development of computer graphics, data visualization and anal- 
ysis have provided both qualitative and quantitative insight about fluid flows [20, 21]. 
Many effective ways of conveying large vector field data into geometric objects, such 
as arrows, glyphs, hedgehogs, streamlines, streamsurfaces and texture mapping, have 
been developed over..the past two. decades. In recent-years, critical point theory based 
topological methods have been introduced into vector field visualization for better 
understanding of the local and global structures of vector fields. Flow visualization 
also serves as the basis for classification of fluid flows [22, 8, 9, 23, 10, 24, 25]. 

Despite all the great developments in vector field visualization, to our knowledge, 
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there hasn’t been any research on quantitative measure of closeness between vector 
fields. The usefulness of this measurement can be seen when comparing computational 
and experimental flow fields under the same conditions for design purposes or for 
the verification of theories. Furthermore, its applicability can be extended to more 
cumbersome tasks such as navigating through a large database while searching for 
similar topology. It might also be used to better quantify the changes in time varying 
data sets. 

1.3.2 Overview of tensor field visualization 

As with vector fields, there are numerous second-order tensor quantities that are of 
interest for analysis and visualization. In fluid mechanics, the Reynolds stress and 
the strain-rate (deformation) tensors are tensors of considerable interest in turbulence 
studies. The alignment of an medium eigenvector 2 of a deformation tensor with the 
vorticity vector and scalar gradients [26] [27] provides a way of vortex detection. In 
solid mechanics, the study of the Boussinesq tensor tells us the force distribution 
inside a solid body. And in general relativity, the basic field equations are described 
by the Einstein tensor and energy-momentum tensor. 

Second-order tensor data sets are at the heart of many engineering and physical 
disciplines. Yet very few methods have been developed to understand and visual- 
ize such data sets due to their complex nature and their large size. In traditional 
approaches, data are displayed using simple local icons depicting, for example, eigen- 
vectors and eigenvalues at discrete locations. Alternatively, scalar components are 
analyzed and ^displayed in ,two-_ or .three-dimensional-form -as an aid- to understanding 
these data. These classical methods are far from effective, but are almost universally 
used to try to extract important information about underlying physical processes and 
engineering principles. More significantly, through this data presentation process, 

2 see sec: hyper 
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valuable information such as continuity of the data and global structure is discarded. 
These difficulties provided us with a substantial motivation to find efficient meth- 
ods for analyzing and displaying them. Thierry Dermarcelle introduced topological 
methods for analyzing and visualizing second-order tensor data [11]. This approach 
preserves the continuity of the tensor field and reveals the key properties of the field 
without any information loss. A new visual representation-hyperstreamline [12] has 
also been designed in order to visualize the topological skeletons of the tensor field. 
However, his topological study is mainly applicable to two-dimensional space and 
concentrates only on simple, isolated degenerate point theory in the field. We live in 
a three-dimensional world and there are also a lot of cases where degenerate points 
are continuous. Therefore, it is important for us to understand the underlying physics 
of this world. 

1.3.3 Objectives 

Our main objectives in this research are as follows: 

• Develop a quantitative measure for feature comparisons of vector fields, 

• Design a good preprocess procedure to simplify the tensor analysis, 

• Study the topological structure of three-dimensional tensor fields, and 

• Find the representations of continuous degenerate points. 

Vector Field Comparisons 

The concept of critical points in vector fields has been well studied and has proven 
to play a key role in vector field topology. However, there hasn’t been any study on 
quantitative measurement of similarities and differences between vector fields using 
topology. In fluid dynamics, simplifications are necessary in the numerical simulation 


10 


process in order to compute analytical solutions of flows over complex bodies in a 
reasonable amount of time. Researchers formulate various models with the hope 
of capturing the essential features of real flows. For example, Reynolds stress and 
velocity fields in turbulent flows are almost impossible to solve analytically. Many 
turbulence models, both linear and nonlinear, have been suggested to better capture 
real flows. 

For example, in eddy-viscosity models, the average velocity gradient is related to 
velocity 

r\TT OTT A 

-ufUj = VtifoT + -fa) - 3 *>i = !> 2 ,3 (1-7) 

where i/ t =(turbulent velocity scale) x (turbulent length scale) and k is the eddy- 
viscosity. This approach to turbulence closure is extremely attractive from a compu- 
tational point of view, especially in terms of numerical robustness. Second-moment 
closure models (originated by Pope [28]) have been pursued extensively over the past 
decade, with recent efforts being exemplified by the models of Speziale et al‘ [29], Fu 
et al’ [30], Launder and Tselepidakis [31], Craft and Launder [32] and Hanjalic and 
Jakirlic [33]. All these models are derived from the exact Reynolds-stress-transport 
equation and are mathematically elaborate but computationally expensive. Recent 
developments in the construction of nonlinear eddy-viscosity formulation, such as 
The Lien-Leschziner Model [34], The Craft- Sug a- Launder Models [35, 36] and The 
Apsley-Leschziner Model [37], have been made in order to combine the simplicity of 
the eddy-viscosity formulation with the superior fundamental strength and predictive 
properties of second-moment closure. 

However, all these models for a large part have been quite arbitrary. In order to 
test the validity of a certain model, a comparison between the computational and 
experimental results is necessary. However, a metric for the closeness between the 
computational and experimental results remains unsolved. Traditionally, a series of 
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charts for each component of Reynolds stress and the velocity vector have been drawn. 
It is perceptually nonintuitive, as well as difficult to show the intrinsic structural 
differences between the two fields. This dissertation fills the void of not having a 
measurement for feature comparison between vector fields and for second-order tensor 
fields. 

Decomposition of .Tensor Fields As a Preprocessor 

Tensor analysis is a very important yet underdeveloped area in visualization. The 
sheer enormous volume of information and the multivariate nature have been the 
major hindrances. How to simplify the analysis and extract a small fraction of the 
data sets without loss of information becomes the key issue. Every tensor field can 
be decomposed into a deviator and a spherical part (definitions to follow). The 
spherical part is an isotropic tensor and therefore remains invariant to a coordinate 
system transformation. As a result, there is no particular need to study this part 
of the tensor (here we refer to it as “the isotropic tensor”), it serves as a uniform 
bias. It will be shown that the deviator of a tensor is parallel to the tensor itself. 
Therefore, their respective eigenvector fields are identical. Furthermore, the locations 
of the respective degenerate points are also identical. This, in turn, means that the 
topology of a tensor field is identical to the topology of its deviator. A deviator 
determines the properties of a tensor field. 

In this dissertation, extraction of a deviator from a tensor field is described as a 
preprocessor for tensor analysis. 

Topology of Three-Dimensional Tensor Fields 

The difficulty with representing tensor data arises from their complexity. At a point in 
space one frequently has to deal with nine or more components of a physical variable 
or a derived quantity such as vorticity, or stress; and there are usually millions of 
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data points available from computation or experiment. In the past, therefore, with 
few techniques available for tensor data analysis and display (often only on a local 
basis), very little information was obtained from complex calculations concerning fluid 
flows, electromagnetic radiation, or field quantities in general relativity. Much more 
effort is usually expanded on the generation of data than on the understanding that 
can be attained from them. In certain fields of physics, this problem is particularly 
acute. For example, in non-Newtonian fluid mechanics (such as the flow of plastics) 
the constituent equations are often not known, and the tensor data describing stress 
and strain are dependent on the model used for these equations. It would be very 
significant to have methods that can analyze the results of these calculations while 
facilitating insight into the underlying physical phenomena; for example, by efficiently 
analyzing and displaying tensor data sets it may be possible to investigate the effect 
of various forms of the constituent equations on the resulting stress and strain tensors 
in non-Newtonian flows. 

Topological representations have proven to be very useful in fluid mechanics as 
well as in engineering, physics and mathematics. The usefulness is derived from 
the simplicity of the representation, while preserving the richness of the data. A 
skillful designer can make very good use of topological skeletons by having an intuitive 
understanding as to how these skeletons relate to physical flow features which give 
rise to large scale effects such as flow separation, or drag and lift on lifting bodies. 

Tensor field topology has not been studied in much detail for the purpose of visu- 
alization. Therefore we have devised methods, based on mathematical foundations, 
that allow the. top.ology.of-Symmetric tensor-fields to-be determined. In two dimen- 
sions, this process is straight forward, and reveals simple geometric pictures that 
represent tensor data locally as well as globally in a very simple manner, usually on 
a surface [11, 12, 13]. In three dimensions, this process is significantly more diffi- 
cult to carry out and theoretical notions are not easily derived. Yet, the results are 


Introduction 


13 


strikingly simple and effective in those cases that we have studied so far. In fact, 
tensor field topology is often simpler than vector field topology, a counter-intuitive 
result given that tensor data appear more complex and rich than vector data sets. 
Three-dimensional tensor field topology will be discussed in detail in Chapter 4. 

Representation of Continuous Degenerate Points 

Previous studies of tensor field topology only deal with isolated degenerate points in 
two-dimensional case [38]. However, points of double degeneracy in three-dimensional 
space may appear along lines or surfaces. In this dissertation, the concept of a 
’’control function” of a deviator is proposed. This function determines the existence 
of degenerate points and assists in defining the lines and surfaces along which the 
degenerate points lie. Following the study of the tensor field in the neighborhood 
of degenerate points, one can display a full representation of the three-dimensional 
tensor field. 


1.4 Visualization Techniques for Vector and Ten- 
sor Fields 

In this section, we discuss a few new developments in visualization techniques and 
show their advantages over the conventional methods. For vector field visualiza- 
tion, the most straight forward method is drawing point icons; while for second-order 
tensor .fields, -the .traditional .method, is-drawing .contour- plots for individual compo- 
nents or glyphs for eigenvectors at discrete locations. All these techniques have their 
limitations for visualizing continuous multivariate data sets. Streamlines (including 
streamtubes and streamsurfaces) and hyperstreamlines are new visual representations 
for vector and tensor fields, respectively. These techniques are designed especially for 
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displaying multidimensional multivariate information while maintaining continuity. 
Texture mapping is a very good way of rendering two-dimensional or 2D slices of 
three-dimensional continuous data sets, both for vector and tensor fields. 

1.4.1 Vector Fields 

Point Icons 

Point icons are useful in visualizing 2D slices of 3D vector fields. There are sev- 
eral candidates such as arrows, wedges and hedgehogs. A comparative study shows 
that arrows are most efficient at conveying vector information from volumetric data 
sets [39]. This method is simple and straight forward, but is impractical when ap- 
plied to the entire volume because of the visual clutter. Therefore, the density of 
displayed arrows must be kept low. However, this make it difficult to comprehend the 
underlying structure of the vector field by mentally interpolating adjacent arrows. 

Figure 1.1 shows the velocity field of the flow past a circular cylinder. Arrows 
represent the velocity vector at various locations. Color maps kinetic energy density. 
Although valuable information is depicted, the display of merely 4% of the arrows 
renders the image overly cluttered and does not convey the intrinsic continuity of the 
data. 

Streamlines, Streamtubes and Streamsurfaces 

Line icons are more efficient in the sense that they provide a continuous representation 
of the data. Consider a vector field u(x, t) at time to, where streamlines are integral 
curves satisfying: 

Ht 

— = v(x,t 0 ) (1.8) 

where x is the position in space and s is a parameter measuring distance along the 
path. Streamlines are everywhere tangent to the flow v(x, to). 
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Figure 1.9: Ellipses with principal axes representing the two eigenvectors £) of an 
elastic stress-tensor field. 

of the axes of the ellipses/ellipsoids represent the two/three-dimensional eigenvector 
fields (eigenvectors and eigenvalues) vl (Figure 1.9) of the tensor field. The advan- 
tage of this representation is that it remains invariant under rotation. However, el- 
lipses/ellipsoids suffer from the the same limitation, namely, visual clutter (Figure 1.9. 
as arrows in vector fields. 

Hyperstreamlines 

The topology of a tensor field T(x) is the topology of its eigenvector fields vfix), 
i = 1.2.3 (16). Line icons representing eigenvector fields were introduced by Thierry 
Delmarcelle [11. 12] and are able to embed the multivariate information of tensor 
fields along trajectories in 3 D space. 

Definition 3 (Hvperstreamline) .4 geometric primitive of finite size sweeps along 
the ..longitudinal eigenvector^ field. V[. while, stretching rin- the -transverse plane under 
the combined action of the two transverse eigenvectors. v t . and v t2 ■ Hyperstreamlines 
are surfaces that envelop the stretched primitives along the trajectories. We color 
hyperstreamlines by means of a user-defined function of the three eigenvalues , usually 
the amplitude of the longitudinal eigenvalue. [38] 


26 



Figure 1.10: 2-D and 3-D symmetric tensor fields. Orthogonal eigenvector fields Vi 
are represented as bidirectional arrows. 

Color and trajectory of hyperstreamlines represent the longitudinal eigenvector 
field and cross-sections encode the transverse eigenvectors (Figure 1.11). Thus, hy- 
perstreamlines fully represent tensor data along continuous trajectories. Similar to 
streamlines, hyperstreamlines are generated by integrating the eigenvector fields. Hy- 
perstreamlines are referred as “major,” “medium,” or “minor,” depending on the 
corresponding longitudinal eigenvector field that defines their trajectories. 

1.4.3 Texture Mappings 

Textures have been well known in flow visualization. The simulation of particle con- 
vection leads to texture [43, 44, 45]. The clouds, smoke and other typical textures are 
perceived in experimental visualization when many particles are used. However, tex- 
ture mapping is a.recently. developed technique-in -the area of-scientific visualization. 
It renders streamlines/hyperstreamlines without integrating the vector /eigen vector 
field and gives a continuous, space-filling view of a 2D field or 2D slices in a 3D field. 
The drawback of texture mapping is that it’s viewpoint dependent. 

Many texture mapping techniques have been proposed. Van Wijk’s spot noise 
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Chapter 4 discusses the topology of 3D second-order tensor fields. First, a de- 
composition procedure is introduced. It serves as a preprocessor to simplify the tenor 
analysis. Second, the existence conditions of trajectories are discussed in Q — R plane. 
Third, the representations of continuous degenerate points in 3D are defined by con- 
trol functions. Finally, the basic structures about simple degenerate points in 3D 
space are discussed. 

Vector and tensor field visualization is a relative new yet very important area. Vi- 
sual representations are most efficient ways for scientists and engineers understanding 
complex physical phenomena. Chapter 5 summarizes the contributions of the current 
research and discusses the future research issues. 


Chapter 2 


Review of Vector and Tensor Field 
Topology 


Topology is the basis for the visualization of vector and second-order tensor fields. The 
purpose of visualization is to help scientists and engineers understand the behavior 
of trajectories in a certain domain. Topology methods seek the surfaces in 3 D or 
the curves in 2D which divide the domain into regions where trajectories behave 
differently. The set of dividing curves or surfaces defines the vector/second-order 
tensor field and is determined by the critical/degenerate points in the field. 

This chapter. isJntended.to give an overview, of .the-state of -the art of the research 
on vector and second-order tensor field topology. It is divided into two parts. The 
first part discusses the critical point theory and vector field topology in both 2D 
and 3D space; the second part discusses degenerate point theory and research on the 
development of tensor field topology. 
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Figure 2.1: An example of a 2D vector field 


2.1 Vector Field Topology 


A vector field V on a subset D is a function assigning to each point P of D a vector 
with its tail at P [17]. Figure 2.1 shows an example in a 2D plane. Placement of 
the vector V(P) with its tail at P is mainly a dramatic visual aid, which adds little 
quantitative information about the field. The essential qualities of the vector V(P) 
are its length and direction. Place all vectors with their tails at the origin, then V(P) 
is described as: 


V(P) = (F(x,y),G(x,v)) 


( 2 . 1 ) 


where F and G are real-valued functions of P = ( x,y ). 
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2.1.1 Critical Point Theory 


A vector field 2.1 determines a system of differential equations: 

i = 

. £ = G(x, 

in some region D. The solutions form a family of directed paths, called integral paths 
of the system, which are tangent to the vector field at each point P. There is exactly 
one integral path that passes through each point P at which V(P) is not zero. The 
picture formed by these paths is called the phase portrait of the system of differential 
equations. The phase portrait is determined by particular points P, called critical 
points 1 , where V(P) = 0 and around which the integral paths gather [17]. 

The most important characteristics of the phase portrait are the number and 
arrangement of the critical points, the pattern of the integral paths about each point, 
and the stability of the critical points. The number and nature of critical points 
won’t change under continuous transformations. These are topological properties 
of the system of differential equations. To characterize the topology, a number of 
definitions and theorems are introduced in the following sections. 



Winding Number 

Consider a continuous vector field V on a closed path L with no zero on L. Starting 
at a fixed point P on L, the vector V(P) will wiggle about during the trip around L 
and return to the point P with some whole number of revolutions (Figure 2.2). 

Proof: The vector V(P) starts at P and after traversing the path L , it returns 
to the same position and with the same direction and length as the original V(P). 
Therefore, the number of revolutions during the trip has to be integers. 

This completes the proof. O 

: see Chapter 1, section 1.2.1 
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Figure 2.2: Winding number of a vector field on a close path L 

Counting these revolutions positively if they are counterclockwise or negatively if 
they are clockwise, the resulting algebraic sum of the number of revolutions is called 
the winding number of V on L , denoted as W (L) [17]. 

Isolated Critical Points and Poincare Index 

Definition 4 (Isolated Critical Points) If for a critical point P, there exists a 
neighborhood of that vector field V in which the vector field vanishes only at P, then 
P is an isolated critical point. 

The Poincare index of V at P, denoted I(P), is defined as the winding number 
W( 7 ) of V on 7 . 

Figure 2.3 shows some examples of Poincare indices of different types of critical 
points. 

. Theorem-1 (The Poincare-Index -Theorem) Let -V be a continuous vector field. 
Let D be a cell and 7 its boundary. Supposing that V is not zero on 7 , then 

W{ 7 ) = /(Pi) + J(P 2 ) + • • • + I(Pn) (2.3) 


where Pi, P 2 , .... P n are the critical points of V inside D. 
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Figure 2.3: Some examples of Poincare indices of different types of critical points. 


The Poincare index measures a property of the vector field that depends only on 
the behavior of V in an arbitrarily small area around P. It is only a crude means of 
classifying critical points, since points of widely different types may have the same 
index (Figure 2.3) [17]. 

Basic Building Blocks 

All critical points are made up of sectors containing one or more of the three types 
shown in Figure 2.4: an elliptic sector, where all paths begin and end at the critical 
point; a parabolic, sector,, where just .one .end -of-the. path is -at the critical point; 
a hyperbolic sector, where the paths do not reach, just sweep past the critical 
point. These three sectors are the basic building blocks of the vector fields surround- 
ing the critical points. The paths that divide each sector from the next are called 
separatrices. 
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Efcrabolic Sector 


Hyperbolic Sector 


Figure 2.4: Three basic sectors 



FniipHr flartry 


Figure 2.5: A typical isolated critical point 
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A typical critical point may have sectors of all three types (Figure 2.5). A critical 
point may have only one type of sector. Those with only parabolic sectors are called 
nodes. Those with only elliptic sectors are called roses, for example, a dipole. Those 
with only hyperbolic sectors are called cross points. Saddle points are cross points 
with four sectors. Still more complicated types are possible with an infinite number 
of sectors, and may be nonisolated critical points [17]. Here, only isolated critical 
points with a finite number of sectors are under consideration. 


2D Critical Points 


Since a vector field vanishes at a critical point, the behavior of nearby streamlines is 
determined by the first order partial derivatives of the vector field. More precisely, 
the 2-D vector field v = (ui,u 2 ) near a critical point is given, in most cases, by the 
first-order expansion 


vi(dx u dx 2 ) ~ %£dx 1 + %£dx 2 
v 2 (dxi,dx 2 ) ^ ^dx 1 + ^dx 2 

where dx 1 and dx 2 are small distance increments from the critical point position. 
Thus, the nearby flow pattern is completely determined by the 2x2 Jacobian matrix 
J, whose elements 


T _ dVi 
13 ~ dXj 


(2.4) 


(i,j = 1,2), are evaluated at the critical point position. 

Different patterns arise that are characterized by the eigenvalues of the matrix 
J. Figure 2,6.showa.how the ^eigenvalues of .X classify- ^critical point as an attract- 
ing/repelling node, an attracting/repelling focus, a center, or a saddle. Real eigen- 
vectors of J are tangent to the streamlines ending at the critical point. A positive 
or negative real part of an eigenvalue indicates an attracting or repelling nature, 
respectively. The imaginary part denotes circulation about the point. 
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1 

T 


Repelling Focus 
Rj and R2 > 0 
I] and I2 0 0 


Saddle Point 

Rj . R 2 < 0 
I] and I2 = 0 


Repelling Node 

Rj and R2 > 0 
Ij and I2 =0 



Attracting Focus 


R ] and R2 < 0 
I] and I2 <> 0 


Center 

Rj and R2 = 0 
lj and I2 < > 0 


Attracting Node 
Rj and R2 < 0 
Ij and I2 =0 


Figure 2.6: 2-D critical points. R1 and R2 denote the real parts of the eigenvalues of 
J, II and 12 the imaginary parts. 

In 2D flows, critical points near the surface of a body known as ’’attachments” 
and ’’detachments” play a similar role as saddle points Figure 2.8. In 2D flows 
near the surface of a body where the velocity is constrained to be zero, the flow is 
mainly tangential and streamlines propagate parallel to the surface. The ’’detach- 
ment” /’’attachment” appear at point where the tangential velocity vanishes and a 
streamline suddenly originates/terminates [53]. 

3D Critical Points 

Simple critical, points of a 2D vector field give rise to a basic. set of patterns shown in 
Figure 2.6. In 3D, critical points can be classified by simple generalization [54] of the 
2D classification shown in Figure 2.6. Alternatively, the classification can be done by 
examining the invariants of the matrix J [22]. 

For a 3D vector field defined over a 3D domain, v = (vi,V 2 ,V 3 ), the Jacobian 
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Figure 2.7: 3-D saddle/saddle/node. 


matrix J is a 3 x 3 matrix whose elements are given by Equation 2.4, for i,j = 1, 2, 3. 
However, J has three eigenvalues and three eigenvectors. Complex eigenvalues always 
occur in conjugate pairs together with a third real eigenvalue. Again, real eigenvectors 
are tangent to streamlines ending at the critical point, and complex eigenvalues denote 
circulation. 

Possible 3D patterns include purely attracting/repelling nodes with eigenvalues 
being all real and all negative/positive, and they behave as 2-D attracting/repelling 
nodes in each of the three planes spanned by pairs of eigenvectors; saddle-saddle-nodes 
with eigen values, heing all real Jbul one.having- a. different-sign, -and they- behave as 
2-D saddles in two planes and as a 2-D node in the third plane; and spiral nodes with 
one real and two complex conjugate eigenvalues, and have an attractive or repelling 
third direction. Figure 2.7 shows an example of an saddle-saddle-node. Algorithms 
for locating and extracting critical points can be found in [55, 56]. 
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Unsteady Flow Fields. Critical points in unsteady vector fields move and even- 
tually merge or split. These phenomena are studied by tracking and representing 
their trajectories over time [57], or by locating interactively critical points in nearby 
space-time regions [58]. 

2.1.2 Global Topology 

Vector field topology depicts schematically the behavior of a large collection of stream- 
lines. The basic constituents of vector field topology are 

• the critical points (Section 2.1.1) and 

• the set of their connecting streamlines (separatrices). 

The topological skeletons divide the field into regions topologically equivalent to 
uniform flow. The representation is highly effective, due to the ease of inferring the 
behavior of every streamline in the flow from these simplified graphs. Also, comparing 
the flow at different time steps is greatly facilitated and can be automated using 
syntactic pattern recognition [20]. 

2D Time-Dependent Fields 

Streamsurfaces described in Section 1.4.1 can be used to visualize two-dimensional 
unsteady flows that depend on time, Reynolds number or any other parameter with 
the third dimension corresponding to the parameter value. Generated by stacking 
instantaneous topological skeletons, streamsurfaces are able to display the topological 
transitions that may occur between consecutive steps. 

Figure 2.9 shows the surfaces in a simulated flow about a 2D cylinder. Time 
increases from back to front. The surfaces are obtained by connecting the streamlines 
and the third dimension corresponding to time. Yellow and blue surfaces correspond 
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Figure 2.8: Topological representation of the 2-D flow past a circular cylinder at two 
different time steps. The flow is coming from the left; at — attachment point; de — 
detachment point; sp = saddle point; ce = center. 

Image courtesy of J. Hetman (Reference [10]). 


to incoming and outgoing separatrices, respectively. Surfaces from attachment points 
are colored orange and those from detachment are colored purple. The periodic 
vortex shedding can be seen in the repeated development and movement downstream 
of saddle-center pairs. 

3 D Separated Flows 

3-D separated flows play a significant role in aerodynamics because separation sur- 
faces are often associated with vortices and recirculation zones, which are important 
structures of the flow far from the body [53]. 

As in the 2-JD_example . of -Figures 2.8 and .2.9, -the. fluid in 3-D separated flows 
moves parallel to the body and then suddenly detaches from the surface, creating 
vortices in the wake. The fluid can also reattach, causing recirculation regions similar 
to the bubble in Figure 2.8. However, detachment and attachment in 3-D flows do 
not arise at isolated points on the surface of the body, but are distributed along entire 
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2D tensor field topology and covered some work in 3D space. We outline some of his 
findings in this section. 

A second-order tensor field T(x) defined on an open subset E of R n is a mapping 
T that associates with each vector x of E a linear transformation onto itself. In 
Cartesian coordinates, T can be represented by n 2 components = 1 ,... , n. 

And the transformation of its components follow the law: 

n,= Y. (2.5) 

P,Q = 1 

under an orthonormal transformation {3 = of the coordinate systems (see Chap- 
ter 1, section 1.2.1). 

The equivalent representation of a tensor field T is using its eigenvectors and 
eigenvalues: 

Vi - A iFi (2.6) 

where i = 1,2, ...,n (n is the dimension) and the eigenvalues are numbered as A x > 
A 2 .... > A n . Eigenvectors are bi-directional due to their sign indeterminacy (Fig- 
ure 1.10). When the tensor components of T are continuous and differentiable func- 
tions of E, its eigenvectors and eigenvalues are continuous and smooth at most points 
(except at degenerate points). 

Any tensor T can be decomposed into the sum of a symmetric tensor S and an 
antisymmetric tensor A: 

S = i(T+.T’) 

A = |(T-T’) 

where T’ is the transpose of T. The eigenvalues and eigenvectors of S are real and 
orthogonal. The antisymmetric tensor has only three independent components that 
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form a tensor known as : ‘axial vector”. 3 The focus in the following discussion will be 
on 2 D symmetric tensor fields with some partial extension in 3 D space. 


2.2.1 Degenerate Point Theory 

Second-order tensor fields have many properties analogous to vector fields. Degener- 
ate points are points where at least two eigenvalues are equal to each other 4 . They are 
basic singularities underlying the topology of tensor fields and play a similar role as 
critical points in vector fields. Away from degenerate points, the hyperstreamlines 5 
are smooth, continuous and non-intersecting trajectories. The tensor data in this 
region are a diffeomorphism to a constant field (Figure 2.12). Therefore, degenerate 
points and their vicinities are the only places that are topologically interesting and 
they determine the behavior of the whole field. 


Tensor Index 

Similar to vector fields, the local behavior about degenerate points can be character- 
ized by a tensor index. The notion of a ’’tensor index of a curve” is first introduced. 
And the concept of ’’tensor index at a degenerate point” is derived from this notion 
which is an extension from vector field topology. 

The ’’curve” under consideration is a Jordan curve: 

3 The antisymmetric tensor field 


is equivalent to the vector field 



in the sense that wi, wa, and uj-j transform as the components of a vector field. 
4 see Chapter 1, section 1.2.1 
5 definition see Chapter 1, section 1.4.2 
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Figure 2.12: Away from degenerate points, continuous and symmetric tensor fields 
are diffeomorphic to constant data. 


Definition 5 (Jordan Curve) A Jordan curve is a curve which is homeomorphic 
to a circle — i.e., a piecewise- smooth, simple, closed curve. 

And the "tensor index of a curve” is defined as below: 


Definition 6 (Tensor Index of a Jordan Curve) Let E be an open set of R 2 and 

let L C E be a Jordan curve. The index /y(Z,) of L relative to a tensor field T £ 
C l (E), where T has no degenerate points on L, is the number of counterclockwise 
revolutions made by the eigenvectors of T when traveling once in a counterclockwise 
direction along L. 


The index /'p(T) can be computed as follows: 

The angle a between the eigenvectors of T and the x-axis (Figure 2.13) is given 


by 


which implies 


tan 2a — 


1 


2T 


■ n 


12 

— T22 

2 Ti 2 


a — — arctan ■ _ 

2 Tn - T 22 


(2.7) 


Equation 2.7 defines a only locally and within an integral multiple of The differ- 
ential of a , however, is a well-defined, smooth, differential form in the whole domain 


E: 



T 22 )(iTi2 — T 12 d(T n — T 22 ) 
(T n — T 22 ) 2 + 4T 1 2 2 


(2.8) 
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Figure 2.13: The two orthogonal eigenvector fields represented as bidirectional 
arrows. 

From Definition 6, the index Tp(T) of a Jordan curve L is equal to the line integral 

1 T (i) = 5?/ i+ da 

where L + means that the curve L is traversed in a counterclockwise direction. Thus, 

J_ f (Tii — T 22 )dT l2 — Ti 2 d(T u — T 22 ) 

2t r/U (Th-T 22 )2 + 4Ti 2 2 

Equation 2.9 will prove extremely useful in analyzing degenerate points (Section 2.2.1). 

Now a tensor index can be associated with degenerate points by enclosing them 
with Jordan curves. 

Definition 7 (Tensor Index at a Degenerate Point) Suppose that T € C l (E ) 
where E is an open subset of R 2 . Let xo G E be an isolated degenerate point of T and 
let L C E be a Jordan curve encompassing Sq in its interior and no other degenerate 
points of T. Then the index of T at the degenerate point x 0 is defined as 

J't’(To) = It( L ) 



The following corollary results from Definition 7: 
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Corollary 1 The tensor index at a degenerate point of a tensor field T € C l (E), 
where E is an open subset of R 2 , is an integer or half -integer quantitrf. 

And there also exists a theorem similar to Poincare index theorem in vector fields 
that deals with a Jordan curve encompassing several isolated degenerate points. 

Theorem 2 Suppose that T € C l (E), where E is an open subset of R 2 containing a 
Jordan curve L. Assume that there are no degenerate points on L but that there are 
a finite number of degenerate points, X\, . . . in the interior of L. It then follows 
that 

I T (L) = yt T (x i ) (2.10) 

J=1 

Basic Building Blocks 

Hyperstreamlines in the vicinity of degenerate points also form several basic building 
blocks as streamlines around critical points in vector field. 

Consider an isolated degenerate point Xq of a tensor field T(x). In most cases, 
the eigenvector fields in the vicinity of xq can be described in terms of three types of 
angular sectors (Figure 2.14): 

• hyperbolic sectors Qj, where trajectories sweep past the degenerate point; 

• parabolic sectors (3j , where trajectories lead away or towards the degenerate 
point; and 

• elliptic sectors y k , where trajectories begin and end at the degenerate point. 

The number of hyperbolic, parabolic, and elliptic sectors are denoted as n k , n p , and 
n e , respectively. Eigenvectors rotate in different directions, when passing through 
sectors of different types along a closed curve that encompasses a degenerate point. 


6 Proof see Reference [38] 
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point. 
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Sector 

Eigenvector Rotation 

parabolic 

counterclockwise 

elliptic 

counterclockwise 

hyperbolic: 

ai > 180° 
< 180° 

counterclockwise 

clockwise 


Table 2.1: Rotation of the eigenvectors during a counterclockwise motion around a 
degenerate point. 

Table 2.1 shows the direction in which eigenvectors rotate during a counterclockwise 
motion around a degenerate point. 

Most degenerate points are built from a series of adjacent hyperbolic, parabolic, 
and elliptic sectors. However, some degenerate points only have one type of sector. 
A “trisector point” (Section 3.3.3) has only hyperbolic sector and a “star-node”, has 
only one parabolic sector. Furthermore, some specific points, namely centers and foci, 
have neither types of sectors. In other words, they are made up of n* = 0 hyperbolic, 
rip = 0 parabolic, and n e — 0 elliptic sectors. 

Hyperstreamline trajectories that separate adjacent sectors at a degenerate point 
play an important topological role. Then “Separatrices” of a tensor field can be 
defined by analogy with a similar concept from vector field topology (References [62, 
63]). 

Definition 8 (Separatrices) A separatrix of a tensor field T is a hyperstreamline 
trajectory which lies on the boundary of a hyperbolic sector at a degenerate point of 

T. 


Separatrices lie on hyperbolic/hyperbolic, hyperbolic/parabolic, and hyperbolic/elliptic 
boundaries, but trajectories along parabolic/parabolic, parabolic/elliptic, and el- 
liptic/elliptic boundaries are not to be considered as separatrices. Also, adjacent 
parabolic sectors are always merged into a unique parabolic sector. 
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2D Degenerate Points 


For a 2D symmetric tensor field, 

( T n (x) T 12 ( x) > 

\ T \2 {%) T22 (%) j 

degenerate points satisfy the conditions: 

TliOeo) ^22(^0) = 0 

< 

T 12 (xb) = 0 

The vicinity of a degenerate point can be expressed as: 


Tn 2 T 22 ~ a{x - x 0 ) + b(y - y 0 ) 
T 12 « c(x - x 0 ) + d(y - y 0 ) 


where £q is the location of the degenerate point and 


a = 
c = 


1 dCTii— 7*22) 

2 dx 

dT 12 1 
dx li 0 


l 1 d(Tn — T22) 

U ~ 2 dy 



x 0 


An important quantity for characterizing degenerate points is 


6 = ad — be 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 


where 5 is invariant under rotation. 

For simple degenerate points (6 ^ 0) at xq , the tensor index is defined as: 

/ T 0?o) = ^sign{5) = ±| (2.14) 

It follows that there are only two patterns possible for simple degenerate points 
(Figure 2.15): 

• Trisector point: 6 < 0, Tj'(xo) = and rih — 3, n e — 0 and n p = 0, 

• Wedge point: 5 > 0, /'p(fo) = \ and n* = 1, n e = 0 and n v — 1 (special case 
when two separatrices are equal, then n v — 0) 
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s i 



Wedge 

points 

5 > 0 

I T = 1/2 



Figure 2.15: Simple degenerate points. 6 = ad — be / 0 and I -p = tensor index. 
Trisector (5 < 0) and wedge (6 > 0) points. Trajectories s* are separatrices. 
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Extensions to 3 D Degenerate Points 

In this section, some extensions of the theory of degenerate points to 3 D space will be 
discussed. In this case, there are three real eigenvalues 7 , A*, and three orthonormal 
eigenvectors, e*, i = 1,2,3, at each point of space. 

Vi(x) - \i(x)ei(x) i = 1, 2, 3 

Various types of degenerate points xq exist: 

A double degenerate point: 

Ai(f 0 ) = A 2 (f 0 ), 

A 2 (x 0 ) = A 3 (£ 0 ), 

and a triple degenerate point: 

Ai(x 0 ) = A 2 (x 0 ) = A 3 (f 0 ), (2-17) 

respectively. 8 

Consider double degenerate case where two eigenvalues are identical — for exam- 
ple, Ai(xo) = A 2 (fo) > A 3 (fo). The tensor field is degenerate in the plane orthogonal 
to u 3 (xo) within which locally two-dimensional patterns such as wedge points, trisec- 
tors, and multiple degenerate points can occur. Figure 2.16 represents the tensor field 
around a wedge point and a trisector in the plane orthogonal to u 3 (xo). 

In fact, the patterns in Figure 2.16 should be drawn on a surface normal to the non- 
degenerate eigenvector field, v 3 (x). However, when close enough to x 0 , this surface 
can be approximated with arbitrary precision by its tangent plane P at Xo . (Away 
from f 0 , trajectories depart from P.) 

Let Ai(fo) = A 2 (x 0 ) = A, B° = {/?£}, i,j = 1,2,3, be the rotation matrix that 
brings the third axis of the coordinate system along e 3 (x 0 ) — the unit eigenvector 

7 Here, we only study symmetric tensors. 

8 Recall that Ai(x) > A 2 (x) > A 3 (x) at every point x. 


(2.15) 

(2.16) 
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Figure 2.16: Three-dimensional wedge, point and trisector .around a degenerate point 
where Ax = A2 > A3. 


56 


associated to A 3 (x 0 ) (Figure 2.16). Transforming T by the rotation B° = {/??•} leads 
to new components 

3 


T'i, = £ W,Tp, 

p,q=l 


At x 0 , T’ is diagonalized — i.e., 


T’(xo) = {:?;.(*„)} = 


0 A 0 
^ 0 0 A 3 ^ 


Let x\ — £p =1 (3i P (x p - x°), i — 1, 2, 3, be the new coordinates. x\ = 0 at the location 
of the degenerate point x 0 . x\ and x' 2 run along the plane P while x' 3 goes in the 
direction normal to P. This new frame of reference is shown in Figure 2.16, where 
we translated the origin away from the degenerate point for clarity. The components 
of T’ in the vicinity of x 0 can be expanded to first-order terms as 

T ;, (/) * + £ a\ jk x' k 

k - 1 


with i, j, k — 1,2,3 and a' ijk = -q£-. Tensor patterns can be analyzed in the plane P 
(Figure 2.16) by letting x' z = 0 and considering the 2x2 block 


/ 


T' 

1 li 

T' 

1 12 

T' 
1 12 

T' 

1 22 


\ 


A + ct'mx'i F a' u2 x' 2 a[ 21 x[ + a[ 22 x' 2 

^ 121*^1 F ^ 122*^2 ^ F <^ 221-^1 F <* 222^2 ) 

The important parameter is 6 = aid! — b'd, where a' = |(o: , ul — a 22 1 ), b' = \{a.' ll2 — 
a 222 ), d = a' 121 , and d' = a' 122 . Trisectors, wedge points, and multiple degenerate 
points correspond to.,5 <C 0,^ >.0, and 5 = 0, . respectively. 

The discussion above concerns degenerate points where Ai(xo) = A 2 (xo) > A 3 (xo). 
The case of degenerate points characterized by Ai(xo) > A 2 (xo) = A 3 (xo) is similar. 
However, fully three-dimensional patterns that exist in the vicinity of degenerate 
points where Ai(xo) = A 2 (xq) = A 3 (xq) = A will be analyzed in Chapter 4. 
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2.2.2 Global Topology 

Tensor field topology describes the behavior of a collection of an infinite number 
of hyperstreamlines in the field. The basic constituents of second-order tensor field 
topology are 

• the degenerate points (Section 2.2.1) and 

• the set of their connecting hyperstreamlines (separatrices). 

The technique to extract the topology of tensor fields and to study topological 
transitions is similar to that of vector fields. Tensor fields are visualized by present- 
ing their eigenvector fields individually. Each eigenvector field is represented by a 
topological skeleton obtained by locating degenerate points and integrating the set of 
their connecting separatrices. 

2D Time-Dependent Fields 

In this section, an example of a stress tensor field in a 2D flow past a cylinder is 
illustrated. The purpose of visualizing this field is to study its topological transitions. 
Fluid elements undergo compressive stresses while moving with the flow. Stresses are 
described mathematically by the stress tensor, which combines isotropic pressure and 
anisotropic viscous stresses. 

O — — pdij + UCij + AOy^— 

where i,j = 1, 2, 3, p = pressure, u* and Vj = velocity, and v and A = 

viscosity coefficients. 

Both eigenvalues of the stress tensor are negative and the two orthogonal eigenvec- 
tors, V\ and v -2 (Equation 2.6), are along the least and the most compressive directions, 
respectively. At a degenerate point, the viscous stresses vanish and both eigenvalues 
are equal to the pressure; degenerate points are points of pure pressure. 
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The first step for understanding the structure of the stress-tensor field consists in 
tracking the motion of degenerate points over time. Degenerate points are located at 
each time step. 

Figure 2.17 shows results at a specific instant in time. The dots mark the location 
of simple degenerate points. The top view shows the overall data while the bottom 
view focuses on the region close to the body. The colored background encodes the 
magnitude A 2 of the most compressive force, from very compressive (red), to mildly 
compressive (orange, yellow, green), to little compressive (blue). Wedge points and 
trisectors are represented as black and white dots, respectively. It is shown clearly 
that wedge points are situated in the wake about the cylinder axis while trisector 
points are located off-axis. The images in Figure 2.17 belong to two video clips that 
visualize the motion of degenerate points over time. 

Topological Skeletons 

Topological skeletons are obtained by detecting degenerate points and integrating the 
set of their connecting separatrices. Trisector points in tensor fields play the topo- 
logical role of saddle points in vector fields. As shown in Figure 2.18, they deflect 
adjacent trajectories in any one of their three hyperbolic sectors toward topologi- 
cally distinct regions of the domain. Wedge points possess both a hyperbolic and a 
parabolic sector. They deflect trajectories adjacent in their hyperbolic sector, and 
terminate trajectories impinging on their parabolic sector. The texture represents 
the most compressive eigenvector of the stress tensor (# 2 ). Color encodes as before 
the magnitude of the . compressive force. (A 2 _), from, most.compressive (red) to least 
compressive (blue). The structure of the tensor field is illustrated by superimposing 
the topological skeleton of V 2 - 

The orientation of the eigenvectors at any point in the plane from the topological 
skeleton can be inferred as follows: hvperstreamlines curve their trajectories so as to 
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follow the shape of the separatrices, bending around wedge points. In Figure 2.18 for 
example, wedge points are located around the cylinder’s axis and trisectors quickly 
move off-axis. The separatrices that emanate from the trisector points delineate the 
global structure of the most compressive eigenvector; hvperstreamlines oscillate about 
the cylinder’s axis, presenting a distorted, wave-like pattern with sharp turns at the 
wedge points. 9 


2.3 Chapter Summary 

This chapter is an overview of vector and second-order tensor field topology. It 
presents the state of art in the theory of vector and second-order tensor field topology. 
Vector field topology has been studied for many years dating back to 19th century 
while second-order tensor field topology has not been studied until recently due to its 
complexity. 

Vector fields are characterized by their critical points and the set of separatrices 
connecting these points. In the vicinity of a critical point, the local topological struc- 
ture is determined by the Jocobian matrix and classified primarily by the Poincare 
index. There are six basic patterns of a simple critical point in a vector field:an 
attracting/repelling node, an attracting/repelling focus, a center, or a saddle. By 
integrating the separatrices emanating from critical points, a topological skeleton can 
be constructed to depict the global topology of a vector field. Topology theory in 
3D space can be studied similarly. The critical points in 3D can be classified in its 
three eigenvector planes. Applications to 2D time-dependent flow and 3D separated 
flows show that visualization technique such as streamlines, streamsurfaces combined 
with topology can capture the key properties such as vortices as well as provide a very 
insightful picture of a vector field. Yet, an important subject as how to quantitatively 


9 A video of this animated flow can be requested from Prof. Lambertus Hesselink 
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measure the closeness between different vector fields has not been addressed. This 
study will be covered in Chapter 3. 

Second-order tensor fields have many properties analogous to vector fields. Degen- 
erate points are basic topological elements and play a similar role as critical points in 
vector fields. In the vicinity of a degenerate point, the local structure is determined 
by its invariant S. Based on the Poincare index, the concept of a tensor index is 
developed and classifies the degenerate points. Because of the sign indeterminacy, 
the simple degenerate points in 2D only have two basic patterns: a trisector point 
and a wedge point. The pattern of a double degeneracy in 3D can be studied simi- 
larly after an appropriate rotation in space. A topological skeleton of a tensor field 
can be obtained by connecting degenerate points with their separatrices. It depicts 
the global topology structure of a tensor field. An example of a 2D stress field in a 
time-varying flow is given: the global structure and its transition are shown. Topics 
concerning topology' in 3D space, like the structure of a triple degenerate point, the 
importance of a deviator and 3D representation of non-isolated degenerate points, 
will be discussed in Chapter 4. 


Chapter 3 

Feature comparisons of vector 
fields using Earth Mover’s Distance 


Vector fields have numerous applications in physical sciences. The concept of critical 
points in vector fields has been well studied and has proven to play a key role in 
vector field topology. However, a formal study on the quantitative measurement of 
similarities and differences between vector fields using topology has not been done. 
In fluid dynamics, simplifications are necessary in the numerical simulation process 
in order to compute analytical solutions of flows over complex bodies in a reasonable 
amount of time. Researchers formulate various models with the hope of capturing 
the essential features of real flows. For example, Reynolds stress and velocity fields in 
turbulent flows are almost impossible to solve analytically. Many turbulence models, 
both linear and nonlinear, have been suggested to better capture real flows (Chap- 
ter 1, section E3.3).. However, themodels-for .a large part have been quite arbitrary. 
In order to test the validity of a certain model, a comparison between the compu- 
tational and experimental results are necessary. However, a metric for the closeness 
between the computational and experimental results remains unsolved. Traditionally, 
a series of charts for each component of the Reynolds stress and velocity vector has 
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been drawn. It is perceptually nonintuitive as well as difficult to show the intrinsic 
structural differences between the two fields. Past study in our group has focused on 
the geometric structure of vector and tensor fields [9] [38] [15]. However, a topolog- 
ical property of a vector field is a property possessed by this vector field and all its 
topological equivalents. The dissimilarity of their appearances sometimes fools the 
human visual system. For example, vector fields 15 and 16 in Figure 3.10, are two 
swirling flows that have the same topological structure yet appear quite different to 
a human observer. Therefore, a quantitative measurement for comparison of vector 
fields is essential. 

Earth Mover’s Distance (EMD), first introduced by Yossi Rubner et al. for fast 
retrieval of similar images in a large data base [64] [65], computes the minimal amount 
of work that must be performed to transform one distribution into the other by moving 
“distribution mass” around. In our vector field analysis, the feature distribution is 
defined as the set of critical points and their corresponding a, /3 parameters (to be 
defined later) of a vector field. After extracting the feature distribution, a set of 
EMDs are computed. Multidimensional Scaling (MDS) was originally used for color 
perception and spatial frequency discrimination [66]. Here, with the help MDS, we 
can display the similarities and differences between selected vector fields. 


3.1 Clifford Algebra 

In 1876, Clifford introduced “geometric algebras” [67] which is now known as “Clif- 
ford algebras”. It is a generalization of. Grassman’s exterior, algebra and Hamilton’s 
quaternions, both of which sought to capture the geometric and algebraic properties 
of Euclidean space [68]. 

Roughly speaking, a Clifford algebra is an associative algebra with unit element 
1 into which a given Euclidean or Minkowski space may be embedded, in which the 
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corresponding quadratic form may be expressed as the negative of a square. By using 
Clifford algebra, the geometry of a complex field can be described as a multiplication 
of simple vector fields. 

Definitions 

Definition 9 (Quadratic Form) Let V be a finite- dimensional vector space over 
the scalar field F, where F = R or C. A quadratic form on V is a mapping Q : V — »• F 
such that 

Q(\v) = x 2 Q{v), x €F, u-ev; 

and the associated form 

B(v, w) = |{Q(u) + Q{w) - Q(v — 10)}, v,w e V, 

is bilinear. 

Definition 10 (Quadratic Space) When a quadratic form Q on V exists, the pair 
(V, Q) is said to be a quadratic space. 

Definition 11 (Clifford Algebra) Let (V, Q) be an arbitrary quadratic space, A 
be an associate algebra over F with identity 1 and v: V — » A an F-linear embedding 
of V into A. The pair (A, v) is said to be a Clifford algebra for (V, Q) when 

(i) A is generated as an algebra by v(v): v € V and Ai: A € F, 

(ii) (i/(u)) 2 = —Q(v) 1, all v € V. 

Remark:, .conditional) is -a minimality, restriction on -the "size’ of A and condition 
(ii) ensures that A is an algebra in which there exists a ’square root’ of the quadratic 
form -Q [68]. 

These definitions are basic to construct a ’’geometric algebra” for vector field 
analysis. For a usual 2D vector space V, a four-dimensional algebra G 2 can be 
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constructed with the basis as: 



( 1 0 ^ 


r — ^ 
1 

0 

( 0 l] 


/ 1 0 N 

1 = 


II 

■CO 


II 

£ 

II 

(N 



1° 1 ) 


l 1 °J 

\ 1 0 j 


K o - 1 ) 


The rules of multiplication are: 


= e j l = e j j = 1,2 

(3.1) 

= e) = l j = 1,2 

(3.2) 

i 2 = -l 

(3.3) 

eie2 = C2 e i = i 

(3.4) 


This kind of construction can be extended to any dimension n, and details can be 
seen in References [68, 69]. 


3.2 Vector field representations using Clifford Al- 
gebra 

In [70, 71], Sheuermann et al. introduced Clifford vector fields which represent vector 
fields in the four-dimensional algebra introduced in Section 3.1. The usual vectors 
( x,y ) G R 2 are identified with 

xe\ + ye 2 G E 2 C G 2 

and the complex numbers a + bi G C with 

ol T bi G G 2 


For a real vector field 


v:R 2 -> R 2 
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(x,y) ^ (ui,u 2 ) 

the Clifford vector field is set as 

E :E 2 E 2 

r = xe i + ye 2 i-> Vie x + v 2 e 2 

A Clifford vector field is just a multivector field with values in R 2 C G 2 

v:R 2 ~+R 2 cG 2 

Let 2 = x + iy, z = x — iy be complex numbers in the Clifford algebra, 
means 

1 / 

X = - (2 + 2) 

y = h (z ~ z) 

We get 

v{r) = ui(ar,y)ei + u 2 (x,y)e 2 

= ® i G (2:+i) '5 (z:_2) ) ei 

-iv 2 (| (2 + 2 ) , — (2 - 2 )) ei 
= £ ( 2 , 2 ) ei 

Generally, a linear vector field can easily be shown to be: 

v{r) = E(z, z)e\ 

= (02 + 62 4- c)ei 


(3.5) 
This 

(3.6) 

(3.7) 

(3.8) 

(3.9) 


where a, b, c c C. 
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Let E : C 2 — » C be a complex polynomial so that v = E(z,z)e i- Let F k : C 2 —>■ 
C, k = 1, . . . , n be the irreducible components of £, so that £'( 2 , z) = n£=i -£c ? then 
an arbitrary polynomial vector field with isolated critical points can be expressed as: 

v(r) = E(z , z)ei 

n 

= I! ( a k z + hz + c k )e 1 (3.10) 

k= 1 

where z* is the unique zero of a k z + b k z + c k . 

3.3 a- (3 space analysis 

In Chapter 2, six basic patterns of simple critical points are discussed. They are 
determined by the eigenvalues of their Jacobian matrix. More convenient parameters 
based on eigenvalues will be introduced in the following section and they actually 
lead to the discovery of two more basic patterns. 

3.3.1 Introduction to a-(3 space 

For a linear vector field v = ( az + bz + c)e 1 , let a = ai 4- a 2 i and b = bi + b 2 i. 
Eigenvalues of the Jacobian around its critical point z 0 are 

Xi=h + y/\a\ 2 -b 2 

and 

a 2 = 61 

Define two new parameters: 

a = (3.11) 


S = sign(\af - fci)\/||a ' 2 - bl\ 


(3.12) 


Feature comparisons of vector Helds using Earth Mover’s Distance 


69 





Attracting Focus 
a < 0 


Caiter 
a =0 


0 < 0 


0 <0 



a <0 a < 0 

0=0 0 >o 


Figure 3.1: 8 basic patterns for simple critical points 


Criteria for basic patterns of simple critical points are: 


0 < 0, a = 0 


C 

0 < 0, a > 0 


RF 

0 < 0,a < 0 


AF 

0 = 0, a < 0 


AS 

0 = 0, a > 0 


RS 

0>O, \0\ > \a\ 


Sa 

\0\ < M : 



0 > 0,a > 0 


RN 

0 > 0, a < 0 


AN 


where C stands for center, RN for node, A A for attracting node, RF for repelling 
focus, AF for attracting focus, St for star and Sa for Saddle, see Figure 3.1. Compared 
to Figure 2.6, two more patterns: attracting and repelling stars are added. After 
careful study, it is found that the original six patterns only cover cases for which 
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0 < 0 and 0 > 0 with 0 = 0 left out. The 0 = 0 case is proved to be important in 
order to get a smooth transition between different vector fields (see section 3.3.2) 1 . 

The usefulness of a and 0 parameters is not limited to the classification of critical 
points. They are essential in vector field comparisons (Section 3.4). Therefore, a new 
space a — 0 is introduced with a and 0 being abscissa and ordinate, respectively. It 
is defined as a positive-definite quadratic space with norm 

\a\ = y/a • a 

for all a e a — 0. And the distance of a from b (a and b are two points in a — 0.) is 

|a — 6| = yj(a — b) • (a — b) 

The basic patterns are now points located at (a, 0) in this new space. This is 
a true metric with normal Euclidean distance and vectors in this space obey all the 
rules defined for a regular 2D Euclidean space. 


3.3.2 Circular relationship between eight basic patterns for 
simple critical points 

Lemma 1 Given a dynamic system governed by a vector field (F(x,y),G(x,y)), 

I = 


Tt = G{x ’ v) 


(3.13) 


Let Xq be a simple critical point in the field, and and A 2 be the eigenvalues of the 
Jacobian matrix of x o- 


J = 


/ 


a b 
c d 


1 There are actually a family of logarithmic node corresponding to 0 = 0 case, but the choice of 
y (see section 3.3.2) is rather random, so we just use y = 0 as a representative here 
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where 


_ dFjxjj ) 
u ~ dx 

_ dG(x,y) 
c — dx 


h _ dF(x,y) 

dy 

, _ dG(x,y) 

so’ dy 


Xo 


Xo 


xo 


Then: 1) if Ai and A 2 are real and distinct, i.e. Ai = a + f3 and A 2 = a — (3 (for 
definitions of a and 0 see section 3.3.1), there exists a nonsingular real transformation 
reducing system 3.13 to the canonical form 

ft ~ ax + py 
f tss fix + ay 

2) if Ai = A 2 , i.e a is arbitrary and 0 = 0, there exists a real transformation reducing 
system 3.13 to the canonical form [62] 


dx 

dt 

dy 

dt 


ax 

px + ay 

where p = 0 if b = 0 and c = 0 in system 3.13, and p 0 if b and c do not vanish 
at the same time, in which case p may be any number. 

3) if X\ and A 2 are complex, i.e. Ai = a + i0 and A 2 = a — i0 (a arbitrary, 0^0), 
there exists a real transformation reducing system 3.13 to the canonical form [62] 


dx 

dt 

dfH 

dt 


ax — 0y 
0x + ay 


Proof: Given a dynamic system 3.13, in the neighborhood of a simple critical 
point xq, the vector field can be approximated as a linear field with 


F(x, y) « ax + by 


and 


G(x, y) ss cx + dy 
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Let the transformation matrix be: 


5 = 


Pll Pl2 
P21 P22 


The new Jacobian is: 


J’ = S3S~ l = 


( a' b' ^ 

d d! 


For case 1), 


Because J’ = S3S 


-1 


J’ = 


Ai 0 


0 A. 


2 / 


J’S = S3, then 


A1P11 A1P12 

A 2 P 21 A 2 P 22 j 


( 


apu + cpi2 bpu + dp\2 


y a P 2 \ + CP 22 bp'i\ + dp 2 2 


Equate left-hand side to right-hand side, 


Pii(o — Ai) + P12C — 0 , p\\b + pi 2 (d — Ai) — 0 ( 3 - 14 ) 


and 


P2i(a — A2) + P22 c — 0 , P216 + p22(d — A2) — 0 ( 3 . 15 ) 


If both b and c are zero, then the original system is already in its canonical form: 
otherwise, assume 0 and the nontrivial solution takes the form: 


P 12 _ Ai— a P22 _ A?— q 

P 11 c ’ P21 c 


(3.16) 


Clearly, the determinant of S is not zero and is a nonsingular transformation. Because 
Ai = a -I- /? and A 2 = a — /3, Equation 3.16 can be rewritten as: 


P12 _ a+ 0 -a 
Pu c 


P22 a— 3 — a 

P21 c 


( 3 . 17 ) 
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Choose pxi = P21 = c, and x and y can be transformed into X = p n x + pi 2 y and 
Y = p 7l x + p?zy- Then — \ 1 X and ^ = A 2 Y Now let X = u + v and Y = u — v, 
then 


du _ 1 . dX dY 

dt 2 dt + dt 

= i(AxX + A 2 Y) 

= ^[{a + P)(u + v) + (a - p)(u - v)] 

= ^(2au + 2/3v) 

= au 4- /3v (3.18) 

Similarly, 

dv _ 1 dX dY 

dt 2 dt dt 

= i(A,X-A 2 K) 

= |[(« + /3)(t* + v)-(a- 0)(u - u)) 

- ^(2/5 u + 2 otv) 

= /3u + av (3.19) 


This new system is in the canonical form and is transformed from the original system 
by a matrix 



with determinant = 0c ^ 0. 

Case 3) can be proved similarly except Ai = a+i/3, X 2 = ct—i/3, and let X = u+iv , 


Y = u — iv. 


Case 2) needs some special care. Since Ai = A 2 = a, it is in its canonical form if 
both b and c vanish, and p = 0. Now suppose that b and c are not both zero, say 
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c ^ 0, then from the first equation of ( 3.16), 


Pi 2 _ Ai — a 
Pu c 


and there exists no Pij satisfying equation 3.14 with nonvanishing determinant. Thus 
choose Y = py with p a random real number, then the new system is 


The trajectory here is a logarithmic node with a phase portrait satisfying y — 
^zln(cx) where c is an integration constant. 


Lemma 1 shows that as p goes from — oo to 0 to oo, there are actually a family of 
logarithmic nodes (p = 0 is a special case and the pattern is a star) corresponding to 
the same set of a and 0 with 0 = 0 (Figure 3.2). If p = 0, the original system is in 
its canonical form and the pattern is determined to be a star. If the original system 
is not in its canonical form, then p ^ 0 and is arbitrary, we will choose p = a for 
convenience. Also, the choice of p does not affect the values of oc and 0. 

Theorem 3 The ratio R = & determines the local topological structure of a simple 
critical point. 

Proof: The proof is divided into three parts corresponding to the three cases in 



(3.20) 


with a nonsingular transformation matrix 



This completes the proof. O 


Lemma 1. 
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Figure 3.2: From (a)-(e), p decreases from positive (a), (b) to zero (c) to negative 
(d), (e) 
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Case 1 ) Ai and A 2 are real and distinct, then the system can be reduced to its canonical 
form: 

= ax + gy 
= PX + aY 
The trajectory can be obtained by solving function: 


dX 

, dt 
dY 
k dt 


dX aX + pY 
dY ~ pX + aY 


(3.21) 


For ft^O, divide both the numerator and the denominator on the right hand side in 
equation 3.21 by a, 


dX X + RF 
dY ~ RA + Y 


(3.22) 


For a = 0: 


dX _ pY 
dY ~ PX 
F 
~X 


(3.23) 


Similarly, in case 3) Xi and A 2 are complex, and the system can be reduced to its 
canonical form: 

= aX-pY 
= pX + aY 

The trajectory can be obtained by solving function: 


dX 

dt 

dY 

dt 


dX aX - PY 
dY ~ PX + aY 


(3.24) 


For a 7 ^ 0, divide both the numerator and the denominator on the right hand side in 
equation 3.24 by a, 

dX X - RF 


dY RA' + F 


(3.25) 
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For q — 0: 


dX _ PY 
dY ~ Jx 
Y 
X 


(3.26) 


The operation for these two cases does not change the equation but uses the parameter 


R instead to characterize its phase portrait. 


For case 2) Ai = A 2 , there are two cases: 


a) p = 0, 


dX _aX _X 
dY ~ olY ~ Y 


(3.27) 


b) p ± 0, 


dX 

dY 


aX 


pX + aY 

Because p is arbitrary ( Lemma 1), choose p — a, 


dX 

dY 


aX 


pX + aY 
aX 

aX + aY 
X 


(3.28) 


(3.29) 


X + Y 

In this case, /? = 0 and the value of a does not change the trajectory, and it can be 
characterized as R = 0. 


This completes the proof. O 

Corollary 2 An alternative quantity to determine the local topological structure of a 
simple critical point is the angle 0 = tan _1 (^). 

Proof: From theorem 3, a local topological structure of a simple critical point can 
be characterize by ^ which is the same as tan0 (Figure 3.3). For 0 < 0 < 2i r, there 
is a unique 0 corresponding to every R. 
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Figure 3.3: The point (a, /3) and its equivalent angle © in a — /? space. 
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This completes the proof. O 


The advantage of using 0 instead of R is that 


• when a = 0, 3 = 1 => R = oc, 0 = f and 


• when a = 0, 3 = — 1 => R = — oc, © = y 

Thus, all the values corresponding to various patterns are finite. The value 0 will be 
used for the rest of investigation. 

Based on their a and 3 values, different dynamic systems have their corresponding 
points located in different regions of a-3 space. From theorem 3, a and 3 themselves 
alone do not determine the trajectory of a dynamic system, the key value is their 
ratio. Therefore, a normalization procedure for a and 3, namely, 


a = 


3' = 


y/a 2 +3 2 

JL 


(3.30) 

y/a 2 +0 2 

preserves the dynamic system. For the rest of the discussion, all the o’s and 3' s have 
been normalized. The usefulness of normalization will be seen in calculating EMD 
(Section 3.4). 


Corollary 3 For normalized points in a — 3 space, the unit circle can be divided into 
(Figure 3. f: 

1) a — 1, 3 = 0 and 0 = 0 , a repelling star or a repelling logarithmic node, 

2) ^-<a<l,0<3< y r an d 0 < 0 < ^7 r, a repelling node, 

3) — ^ < a < ^ < 3 < 1 and < 0 < |7r,a saddle, 

4) — 1 < a < — 0 < 3 < ^ and |7r < 0 < n, an attracting node, 

5) a = —l, 3 = 9 and 0 = n, an attracting star or an attracting logarithmic node, 

6) — 1 < ot < 0, — 1 < 3 < 0 and n < 0 < |7r, an attracting focus 

7) a = 0, 3 — ~ 1 and 0 = |7r, a center, 

8) 0 < a < l, — 1 < 3 < 0 and |7r < 0 < 27 r, a repelling focus 
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P 


A 



S 


a 


Figure 3.4: Basic patterns for critical points in a -0 space; £7 for center, RN for node, 
AN for attracting node, RF for repelling focus, AF for attracting focus, St for star, 
Sa for Saddle and R for regular point 


3.3.3 Simple critical points in a-/3 space 

For a vector field with one simple critical point at xq, it can be approximated as 


v = (v X x(x - x 0 ) + v xy (y - y 0 ), v yx (x - x 0 ) + v yy (y - y 0 )) 2 (3.31) 


where 



(3.32) 


Also, in Clifford algebra, this field can be represented by 


v = (az + bz + c)e \ 


where a, 6, c € C and z — x + iy, z = x — iy. Using the property: 


e\ = l, e T e 2 = i = — e 2 e! 


2 see Chapter 2 
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and let o = Oi 4- a 2 i, b = 61 + b 2 i and c — Ci + c 2 i where a\,a 2 , 61 , 62 , Ciandc 2 e i2. 

v = (ax 4- a 2 i)(x 4- iy) 4- 4- M)(z - iy) 4- c x 4- c 2 i 

= (aix - a 2 y + bix 4- b 2 y 4- Ci) 4- (ait/ 4- a 2 x - biy -Pl^x + c 2 )i 
= {ci\x - d 2 y + bix + b 2 y + ci)e? + (ai?/ 4- a 2 x - b x y 4- fox 4- c 2 )e 1 e 2 
= [(fli* - a 2 y 4- b x x 4- b 2 y 4- - (ait/ 4- a 2 x - b Y y 4- b 2 x + c 2 )e 2 ]ei 

Therefore, 


V + ^xj/ 2 / - (fxx^o 4- v xy yo) = (d l 4- b x )x 4- (b 2 - a 2 )y 4- c x 


Vy X x + Vyy - ( Vy X x 0 4" u j, y t/o ) = -(<*2 4" 62 )* 4- (&i - ai)t/ 4* c 2 


Equating all the x and y terms: 


fli 4- &i — v xx 

b\ Gl — Vyy 

(®2 + b 2 ) = 

b 2 G 2 — ^xj/ 


and 


Vi (^xx-^-o 4- v X yyo) 


c 2 {v yx Xo 4- v yy yo) 


Solve equations 3.33, 


2 {.Vxx Vyy) b\ — 2 (.Vxx 4 ” Vyy) 
a 2 2 (Vyx 4“ V X y) b 2 = 2^V X y V yx ) 


(3.33) 


(3.34) 


From equations 3.11 and 3.12, a pair of a, 0 values are readily given. After a 
normalization procedure 3.30, the corresponding point of the vector field 3.31 can be 
located on a unit circle in a- (3 space (Figure 3.4). 
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3.3.4 Regular points in a-(3 space 

If a vector field on E does not have a singular point, then it is diffeomorphic to a 
constant field 

v = const • e\. 

where a = 0 and b = 0. This gives a = 0 and /3 = 0 which corresponds to the origin 
in the a-f3 space. 

3.3.5 Multiple critical points in a-/3 space 

An arbitrary vector field can be represented as a multiplication of linear fields (Sec- 
tion 3.2): 

v(r ) = E(z,z)ei 

n 

= JJ ( a k z + b k z + c k )e i (3.35) 

fc=i 

where z k is the unique zero of a k z + b k z + c k . Around each simple critical point x k , 
a vector field can be also expressed as v(x) = J k (x — x k ) using its Jacobian J*. In 
section 3.3.3, a conversion from a vector field in regular space to a Clifford vector 
field is given. All the a k s, b k s and c k s, therefore, can be computed separately at 
each critical point z k . And a complex field with multiple critical points corresponds 
to multiple points (a x , /3 X ), (a 2 , f3 2 ), . . . , (o: n , /3 n ) in a-/3 space. For example, 

v = z 2 e\ 

is a dipole which has two points both located at (1, 0) in a-/3 space (Figure 3.5). 

3.3.6 The Growth of a Vector Field 


n 

^( r ) = II ( a * Z + b k 2 + °k) e 1 


k=z 1 


A Clifford vector field 
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where n is the total number of critical points in this field. 

This energy here is a quantity that characterizes the critical points of a vector field. 
It is different from the physical energy. The concept ’’work” is used to measure the 
energy differences between two vector fields or the amount of energy used to transform 
one vector field into the other. 

Definition 14 (Work) For two vector fields with feature distributions 


{ (® 1 > P\)i @2)1 - ■ ■ 1 (On; Pn) } 


and 


{( a lJ Pi): ( a 2 i P 2 )) • • • 1 ( a m Pn)}‘ 


The amount of work necessary for transforming one vector field into the other is 
defined as: Work = ^/£"=i((<*i - a-) 2 + (A - P-) 2 ). 


Intuitively, given two feature distributions, one can be seen as a set of discrete point- 
objects with a certain amount of mass of earth spread in space, the other as a collection 
of holes in the same space. The work measures the least amount of energy needed to 
fill the holes with earth and is called the Earth Mover’s Distance (EMD). 

Computing the EMD is based on a solution to the old transportation problem 
from linear optimization. [72] This is a bipartite network flow problem which can be 
formalized as the following linear programming problem: Let I be a set of suppliers, 
J a set of consumers, c*,- the cost to ship a unit of supply from i € / to j £ J 

Cij = y/ fa - oy) 2 + - fij) 2 

and it is the same as the Euclidean distance = ||Ej — Uj|| in a-0 space. A critical 
point either exists as a whole or does not exist, it can not be split. Therefore, the 
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flow fij can only be 0 or 1. We want to seek a set of fa that minimizes the overall 
cost: 

£’MD(x, y) = min^2 Y djhj (3-40) 

id id 

subject to the following constraints: 


fij >0 i € /, j € J (3.41) 


Yj fij — Vii J ^ J 

(3.42) 

id 


^ ! fij — -Zi , i E I 

(3.43) 

jeJ 


YVj = 

(3.44) 

jeJ id 



Where x* is the total supply of supplier i and yj is the total capacity of consumer j. 
Constraint ( 3.42) allows shipping of supplies from a supplier to a consumer and not 
vice versa. Constraint ( 3.43) forces the consumers to fill up all of their capacities 
and constraint ( 3.44) limits the supply that a supplier can send as a total amount. 
Constraint ( 3.44) is a feasibility condition that ensures that the total demand equals 
the total supply, in other words, the distributions have the same overall mass and the 
EMD is a true metric. 

For a set of vector fields, they don’t always have the same number of distributions. 
In order to satisfy constraint ( 3.44), we can create constant fields 1 with a = 0 and 
f3 — 0 to make the supply equal to the demand without changing the vector fields. 
For example, the supplier field is 


Vi = (2z.+ (l.+.i.)z + (— 16d-.8z))((l.d-i)z4-.z-Jr(4 — 4i).)(z4-2iz+ (— A+\2i))e\. (3.45) 

The critical points and their corresponding a, (3 values are listed in Table 3.1. 
There are three critical points in this field (Figure 3.7). And the consumer field is 

i>2 = (iz + (4 + 4 i)){z + 2 z + (12 + 4z))(2z + ( — 8 — 8 i))(z + (4 — 4z)) 2 ei (3.46) 
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Order of critical point 

Location 

i*,P) 

Simple 

(4,4) 

(0.5,0.866) 

Simple 

(-4, -4) 

(0.5774,0.8165) 

Simple 

(-4,4) 

(0, —1.0) 


Table 3.1: Locations of three critical points and their corresponding a, /? values; Order 
of critical point: the order of the lowest non-zero terms in the Taylor expansion of a 
critical point. 


Order of critical point 

Location 

(*,/3) 

Simple 

(4,4) 

(0.0, 1.0) 

Simple 

(-4, -4) 

(0.8944, 0.4472) 

Simple 

(-4,4) 

(1.0, 0.0) 

multiple 
(2nd order) 

(4, -4) 

(0.0, 1.0) 
(0.0, 1.0) 


Table 3.2: Locations of four critical points and their corresponding a, (3 values; Order 
of critical point: the order of the lowest non-zero terms in the Taylor expansion of a 
critical point. 


This field has four critical points (Figure 3.8). Three of them are simple ones; 
however, the 4th one is a 2nd order point which means that there are two identical 
pairs of a , /3 values corresponding to this point. Therefore, there are five points in a-(3 
space associated with these four points. The critical points and their corresponding 
a , (3 values are listed in Table 3.2. The supplier side has three points in the a — (3 
space while the consumer side has five. Now let 


V\ — (2z T (1 T i)z -f- 1)((1 T i)z + z i)(z -I- 2 zz) • 1 • lex (3.4 ( ) 

and the .yector. field . remains .unchanged— However, now we. have two more regular 
points corresponding to 1 with a — 0 and (3 — 0, and both the supplier and the con- 
sumer have 5 points in their feature distributions. All the conditions are satisfied, and 
we are ready to compute the EMD for these two fields and find out the dissimilarity 
between them. 
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3.5 Display of EMDs for a Large Set of Vector 
Fields 

The above discussions are for comparison of a pair of vector fields. If there exist 
a large set of vector fields and we want to compare their topologies, it is necessary 
to display them in a more meaningful way than a ID list sorted by their EMDs. 
Yossi Rubner et al. has used Multidimensional Scaling Method (MDS) [64, 65] to 
display a set of images on a 2D map. Given n objects in a high dimension, the MDS 
method computes a configuration of them in a lower dimension space such that the 
distance between every pair of objects in this low dimension space matches best to 
their real distance in the high dimension. Inspired by their work, w r e compute the 
EMDs between every pair of vector fields and position the vector fields on a 2D map 
such that the distances between the vector fields match their EMD values as close as 
possible. 


3.6 Applications 

Using a 2D map for displaying a set of vector fields can assist navigation in the space 
of vector fields. Computing the EMDs between pairs of selected vector fields and 
positioning them in a 2D map give us a better way to display the query results. 

A simple example is shown in Figure 3.9. A list of simple critical points belonging 
to eight basic patterns are sorted by MDS to show the relationship (distances) between 
them. The.result. in -Figure .3.9. shows .that, all -the critical points -fall into the right 
section in a-0 space and are well matched to Figure 3.4. The sign of 0 controls 
whether the flow is circular or asymptotic. For 0 > 0, the vector field is asymptotic 
and the pattern is either a node or a saddle. For 0 = 0, it forms a star pattern. And 
for 0 < 0, it is circular and the pattern is a focus or a center. The sign of a controls 
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3-6.1 Unsteady Flow Past a 2D Cylinder 

The feature comparison methods can also be used for comparing the topology between 
vector fields in different time steps. For example, Figure 3.13 shows an incompressible 
viscous flow past a 2D circular cylinder 3 . When the Reynolds number is higher 
than 40, the flow is unstable. Vortex shedding occurs behind the cylinder in the 
flow downstream. The alternatively shed vortex pattern is called a Karman vortex 
street [74]. Figure 3.12 shows the computations of EMDs between vector fields in 
different time steps with time step three. The topological patterns of the flow vary 
in a cyclic pattern due to the vortex shedding from the cylinder. The EMD for time 
step 3 and 18 is zero which means that the topology of the flow at these two time 
steps are indentical. This is shown in Figure 3.14. The only difference between these 
two fields is that the flow at time 18 is flipped upside down comparing to that of time 
3. In Figure 3.12, the largest difference in one cycle is between time 18 and time 24 
(or time 3 and time 8) where time 18 only has one center and time 24 has two centers 
and one saddle. This is shown in Figure 3.15. 

3.7 Chapter Summary 

Visualizing vector field topology is a very important subject in vector field visualiza- 
tion and has received much attention. However, to our knowledge, almost no work has 
been done on quantitative measurements for vector field comparisons. In this chap- 
ter, comparisons are based on feature distributions of vector fields. Critical points 
are key features of vector fields. . Simple, critical points are characterized by the eigen- 
values of their Jacobian matrix. A new set of parameters a and (3 are introduced to 
replace the eigenvalues as criteria for critical point classification. The basic patterns 
of simple critical points can be classified into eight categories. A complex field can be 


3 details about the computation of this flow see Reference [73] 
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0.8 





Figure 3.12: The EMDs for flows in different time steps with time step 3. The EMD 
for time step 3 and 18 is zero, which means that flow at these two time steps have 
the same topology. 
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Chapter 4 

Topology of 3D Tensor Fields 


Tensor field visualization is a very important area in multivariate multidimensional 
visualization. Due to their complex nature, tensor data sets usually contain huge 
volumes of information. Although they are seemingly random and chaotic, topology 
study shows that they are nicely structured and intrinsically correlated. Extracting 
the topological skeleton of a tensor field reveals its intrinsic geometric structure and 
captures the nature of the field. The topology of 2D tensor fields has been studied by 
Thierry Delmarcelle [38], this chapter intends to explore further into 3D space. The 
present research to our knowledge is original. 

Degenerate points are the basic constituents of tensor field topology (Chapter 2). 
A thorough study of the behavior of a tensor field in a close neighborhood of its 
degenerate points can lead to a simple topological skeleton that connects these points. 
Because the integral lines of eigenvectors (or hvperstreamlines) in a tensor field never 
cross each. other, except at. degenerate-points,, one can reconstruct the whole tensor 
field based only on a small fraction of the data set. In this manner, we are able to 
compress the data substantially with no information loss, but it should be noted that 
this compression is different from the usual approach used in computer science. Also, 
by displaying only the topological skeleton, we can avoid the visual clutter yet still 
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reveal the essential features of the field. 

Studies on first-order (vector) and second-order tensor fields show that their 
topologies have striking similarities. One wonders why a point with the same eigen- 
values in a tensor field behaves in a similar fashion to a singular point (a point where 
the magnitude vanishes) in a vector field. It will be explained in this Chapter. 

Every tensor field can be decomposed into a deviator and a spherical part (def- 
initions to follow). The spherical part is an isotropic tensor and therefore remains 
invariant to a coordinate system transformation. As a result, there is no particular 
need to study this part of the tensor (here we refer to it as “the isotropic tensor”). It 
will be shown that the deviator of a tensor is parallel to the tensor itself. Therefore, 
their respective eigenvector fields are identical. Furthermore, the locations of the re- 
spective degenerate points are also identical. This, in turn, means that the topology 
of a tensor field is identical to the topology of its deviator. 

In order to study the topology of degenerate points, it is critical to locate these 
points first. In this chapter, a method for locating isolated triple degenerate points 
will be presented. Cases of degenerate points of double degeneracy are more com- 
plicated, they may appear along lines or surfaces. In this chapter, the concept of a 
“control function” of a deviator is proposed. This function determines the existence 
of degenerate points and assists in defining the lines and surfaces along which the 
degenerate points lie. Following the study of the tensor field in the neighborhood of 
degenerate points, one can display a full representation of the 3-D tensor field. 


4.1 Decomposition of a tensor field 

A tensor field in the real physical world is often very complex, therefore reducing it 
into a simpler form becomes very appealing. In this section, the tensor is decomposed 
into a “deviator” and an “isotropic tensor” at each point in the tensor field. 
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Definition 15 (Deviator) A tensor is a deviator D iff it is trace free, i.e., Trace ( D ) 

0. 


Definition 16 (Isotropic Tensor) A tensor is isotropic iff Uij = v5ij, where v is 
a stretch factor. 


Any given tensor T, can be decomposed into 1 : 


where 


and 


Therefore, 


^ T n ( x , y, z ) Ti 2 (x, y, z) T 13 (x, y, z) ^ 


T 2 i {x, y, z) T 22 (x, y, z) T 23 (x, y, z) 
T31 (x,y, z) T 32 (x, y, z) T 33 {x,y,z ) 


^ D n (x, y, z) D l2 (x, y, z) £>13 (x, y, z) ^ 


+ 


£>21 (a:, y, z) D 22 (x, y, z) D 23 (x, y, z) 

^ £>31 (x, y, z ) D 32 (x, y, z) D 33 (x, y, z) ) 
Uu{x,y,z) 0 0 

0 Un(x,y,z) 0 

0 0 £33 (x, y, z) 


D ii ~~ Tii q Tjj 

6 j= 1 


U„ = \t Tji- 

6 3=1 


(4.1) 


Du + D 22 + D 33 = 0 


and 

Du = £22 = D33 

The analysis and examples in the rest of the chapter are all for 3-D space, but the theory holds 
for 2-D space as well 
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Figure 4.1: Eigenvector e 3 is aligned to X3, a local 2D tensor on the plane. 

at any point in the field. Here the deviator and the isotropic tensor are denoted by 
D and by U. 

A tensor can also have a “local 2-D deviator” . Any given tensor can be transformed 
into a space so that at least one eigenvector, for example el, is along one of the axes. 
In the plane perpendicular to el, the tensor is locally 2-D (Figure 4.1). It takes the 
form: 


( r pl 

J n 

T' 

-*12 

0 N 

T' 
1 12 

T' 

J 22 

0 

0 

0 

) 


The locally 2-D tensor can be decomposed into a deviator 


and an isotropic tensor 


D' = 


U' = 


D[ 1 D[ 


0 


12 


D [ 2 Di, 


22 


U[ 1 0 


C /' 2 
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the eigenvalues at these points are equal to zero. This means that these points are 
stress free, a fact that can be verified by an examination of the stress equations. We 
have therefore acquired physical insight into the stress tensor field just by examining 
a basic topological feature, a point of triple degeneracy. Also, after searching through 
the whole field, we found that there are a continuous line of double degenerate points 
inside the solid. The stress tensor has a singular ’’local deviator” along this line. 

Figure 4.2 shows the minor eigenvector field (most compressive force). The ar- 
rows represent the two forces, the two balls indicate the location of singular points. 
The hyperstreamlines on the top surface are the separatrices of the singular points 
which describe the topological structure in their vicinity. Inside the solid, the hyper- 
streamlines form trisectors of double degenerate points which indicate the pattern of 
topological structure in the vicinity of their singular local 2-D deviators. Colors of 
the hyperstreamlines encode the magnitudes of minor eigenvalues. We also varied the 
strength of the two forces . It turns out that the stress free points move along a circle 
as the force ratio is varied. Video clip 1 shows the motion of a double degenerate 
point inside the solid with variation of the two forces; Video clip 2 shows the motion 
of two triple degenerate points on the top surface(Figure 4.2) with variation of two 
forces; Red balls indicates the degenerate points. 


4.2 Physical meaning of a deviator and an isotropic 
tensor 

4.2.1 Isotropic Tensor 

Based on section 4.1, one can see that the tensor U is an identity tensor multiplied 
by Ua- Therefore, at any given location in the field, it behaves the same in every 
direction, in other words, it is an isotropic tensor. Since it is isotropic throughout the 


no 


whole field, it is of no particular physical interests to the following study. 


4.2.2 Deviator 

The deviator, in contrast to the isotropic tensor, has a different behavior in all 3 
principal directions except at a singular point where all of its components are zero. 
It is indeed the deviator that represents the deviations from the originally isotropic 
field and makes the entire tensor field so complex and diverse. It is then reasonable 
to think of a real tensor field as a deviator superimposed onto an isotropic tensor. 
By subtracting the contribution of the isotropic tensor from the tensor field, the 
deviator becomes dominant and it enables us to clearly show the topology and the 
fluctuations of the field without the disturbance of the sometimes dominant isotropic 
contribution. One can think of this as a way to reduce the effect of the constant 
background, while making the variations from constancy more pronounced. The 
major focus is then on the deviator, as all the following sections are dedicated to 
the analysis of the topology of the deviator. 


4.2.3 Stress Tensor vs. Viscous Stress Tensor — An Example 


This section shows an example of different effects of an isotropic tensor and a deviator 
in a flow field. The flow' is incident on the front of the hemisphere cylinder and 
proceeds to the back of the body at a skewed angle. The viscous stress tensor in the 
field is expressed as: 


,dv t 

a '> ~ ^dT, + 

where i,j = 1,2,3. And the stress tensor is: 


dxi’ 


£7jj P$ij T <J{j 

where P is the pressure. The only difference between these tw'o tensors is that the 
stress field has a uniform pressure component. Pressure in a stress tensor comes 
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from outside forces acting on the medium and it can be disturbingly large. In Fig- 
ure 4.3, hyperstreamlines (Chapter 1) are integrated along major eigenvectors, and 
the medium and minor eigenvectors residing in the cross sections. The size of the 
tube is determined by the values of the medium and minor eigenvalues. Because the 
pressure component in this field is rather large, some small variations from the vis- 
cous part won’t make an impact on the appearances of the tubes and they all appear 
circular. This makes it difficult to see the actual effect from the flow itself. However, 
in Figure 4.4, the isotropic contribution is subtracted from the tensor. Therefore, 
pressure no longer has any effect and the resultant tensor is a deviator. The hyper- 
streamlines are now all anisotropic, which means that the viscous stresses are very 
different in the directions of the medium and minor eigenvectors. From the compar- 
ison of these two figures, we can see that a deviator truly reflects the nature of a 
tensor field. 


4.3 Degenerate Points in a Deviator 

Because the trace of a deviator is zero, degenerate points in a deviator take on some 
special features. 

Definition 17 (Singular Point) A singular point in a tensor field is a point where 
all eigenvalues of a tensor vanish, in mathematical representation, it is a zero matrix. 

For a deviator (of a symmetric tensor), 

D\\{xq) Di 2(0:0) Diz{x$) 

^12(^0) D22 (3^0 ) D 2 z{Xq) 

K D\z(xq) D 2 3 ( 2 : 0 ) Dzz{xq) y 
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a triple degenerate point satisfies the conditions: 

Dn — D22 = 0 
D22 — D33 = 0 
' D12 = 0 

D13 = 0 

D23 = 0 

From the first two conditions, we get 

Dn = D 2 2 = D33, 


and also 

Du + D22 + D33 = 0 , 


therefore, 


Dn — 0 D22 — 0 D33 — 0 

So, all the components of a triple degenerate point are zero and this is a singular 
point similar to a critical point in a first-order tensor field (vector field) where all the 
components vanish. A double degenerate point in a local deviator 

' On D' u 0 ' 

D12 D22 0 

^ 0 0 A3 j 

also brings 


D[ 1 — 0 D '22 — 0 D' 12 — 0 

This corresponds to a point where only two components vanish (0, 0, V3) in a vector 
field. However, both of them are dependent on the coordinates. These features of 
triple and double degenerate points shed some light on the similarities of topological 
properties between the first- and second-order tensor fields. 
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4.4 The Existence Conditions of Hyperstreamlines 
in Q-R Space 


The eigenvalues of a 3-D tensor T determine the existence of its trajectories and can 
be obtained by investigating its characteristic equation. 


T(A) 


A -Tii 

-P 12 


-r 12 

A — 

“T 23 

-T 13 

— T 2 3 

a-t 3 


= \ 3 + p\ 2 + q\ + r 


where: 


P — T\\ + T 22 -(- T 33 


(4.3) 

(4.4) 


Tu 

T 12 

+ 

Txx 

T13 

+ 

T 22 

T23 

t 12 

t 22 


T13 

T33 


T 2 3 

T33 


T11 

T\2 

T 13 

Ti 2 

t 22 

T‘2Z 

Tl 3 

T 2 3 

T33 


The coefficients P, Q and R are all tensor invariants. 


(4.5) 


(4.6) 


Theorem 5 In P-Q-R space, hyperstreamlines only exist between surfaces S\ and S 2 , 


which are, respectively, given by: 


2 P 3 +9PQ +.2 (P 2 - 3 Q) 3 ' 2 
27 


+ R — 0 


(4.7) 


and 

2 P 3 + 9 PQ - 2 (P 2 - 3Q) 3/2 


27 


+ R — 0 


(4.8) 
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Figure 4.5: (a) v4(A) is a general cubic polynomial with three real distinct roots; (b) 
^4(A) has one double root and one single root: Ai = A 2 > A 3 (c) .4(A) has one double 
root and one single root: Aj > A 2 = A 3 ; (d) A(X) has one triple root: Ai = A 2 = A 3 . 
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Proof: Since the tensor T is real and symmetric, all eigenvalues must be real. Its 
characteristic function A(A) is a cubic polynomial, therefore, it must have three real 
solutions (Figure 4.5). When A{ A) has double roots, it satisfies the conditions: 






A( A) = A 3 + PA 2 + Q\ + R = 0 
^ = 3A 2 + 2PA + Q = 0 


(4.9) 


Solving equations 4.9 gives two condition equations 

„ 2P 3 + 9PQ + 2 (/* - 3Q) S/2 , „ „ 

E. = 57 + " = ° 


(4.10) 


and 


r. 2P 3 + 9PQ - 2 (P 2 - 3Qf 2 , „ „ 

E» = ^ + « = 0 


(4.11) 


corresponding to Ai = A 2 > A 3 and Ai > A 2 = A 3 , respectively. Therefore, 


• when Ei < 0 and E 2 > 0, there are three real distinct roots (Figure 4.5 (a)), 

• when E\ = 0 and E 2 > 0, one double root and one single root exist (Figure 4.5 

(b)), 

• when Ei > 0 and E 2 = 0, one double root and one single root (Figure 4.5 (c)) 
exist and 


• when Ei = E 2 , one triple root (Figure 4.5 (d)) exist. 

From Figure 4.5, obviously if E x > 0 or E 2 < 0 then A(A) will have only one real root 
and a pair of complex conjugate roots. Therefore, the conditions for all A’s to be real 
or hyperstreamlines to possibly exist are that Fj < 0 and E 2 > 0. In P-Q-R space, 
this is equivalent to saying that hyperstreamlines only exist in between surfaces: 


E = 2P 3 + 9PQ + 2 (P 2 — 3Q) 3 ^ 2 + R = Q 


and 


27 


e 2 = 2f3 + 9p g- 2 ( p2 - 3< ?) 3/2 + R = o 


(4.12) 


(4.13) 
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Figure 4.6: Q"-R‘ plane 


This completes the proof. O 

On surface Si, A a (To) = A 2 ( To ); and on surface S 2 , A2 (To) = A3 (To)- Triple 
degeneracy Ai (To) = A 2 (To) = A 3 (To) occurs at points where S- L and S 2 meet, i.e. 

R=§,Q = : ¥• 

Since the tensor and its deviator have the same topology-, it is sufficient to examine 
only the deviator D. By definition, 

Du — D'y) -r D33 = 0. 


and therefore, P = 0. The coefficient Q can also be presented as: 

O' = -\(Dn + ^1,+ Df 3 + 2 D? 2 + 2D? 3 + 2 £fe) (4.14) 

The characteristic equation now becomes: 

A 3 + Q’X + R = 0 (4.15) 

And the P-Q-R space reduces to the Q‘-R‘ plane. 

Corollary 4 In Q’ — R' space, the hyperstreamlines exist between curves L\ and L 2 , 
which are. respectively, given by: 

R!_ _ f-Q’) 1 


(4.16) 
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and 



3 


W_ 

2 


(4.17) 


Proof: From theorem 5, the hyperstreamlines only exist between surfaces: S\ and 
S 2 . Now substitute P — 0 into equations 4.10 and 4.11, Si and S 2 reduce to curve 
L\ and L 2 on the Q' — R' plane(Figure 4.6), given respectively by 


This completes the proof. O 

Double degeneracy occurs on the curves Li and L 2 and triple degeneracy occurs 
where Li and L 2 meet, i.e. Q' — R' = 0. 

4.5 Control functions of the singular points 

From section 4.4, the characteristic equation of a deviator takes the form: 

A 3 + QX + R = 0. 

The standard form for the roots of this cubic equation is: 



(4.18) 


and 



(4.19) 
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Where 


and 



This equation has three distinct real roots when A < 0 and multiple real roots 
when A = 0. This study deals with symmetric tensors, and therefore the eigenvalues 
are always real. This, in turn, means that for all regular points in the tensor field 
A has to be less than 0. The quantity A = 0 only at a degenerate point, where the 
characteristic equation gives multiple roots. 


• For double roots, (^f ) 2 = — (^) 3 / 0, it is a double degenerate point, 

• For triple roots, Q ~ R — 0, it is a triple degenerate point. 


The existence conditions derived in section 4.4 can also be obtained from A (x, y, z ) = 
0, but from a different point of view. Function A (x, t/, z) controls the occurrences of 
degeneracy of a deviator, thus it is called a “Control Function”. 


4.6 The nature of a degenerate point 

Theorem 6 A tensor is of triple degeneracy iff its deviator is singular. 

Proof: If a tensor T is of triple degeneracy, then T and its isotropic tensor U 
are identical. Therefore, its deviator D = T — U = 0, or D is singular; on the other 
hand, if-D is singular,- then T = D + U = T. 

This completes the proof. O 


Theorem 7 A tensor is of double degeneracy iff its deviator is nonsingular and one 
of its local 2-D deviators is singular. 
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Proof: If a tensor T is of double degeneracy, let its eigenvalues be Ai, A2 and A 3 , 
and also Ai = A2, then it has one distinct eigen direction e$ associated with A3. We 
rotate the coordinate system so that one axis is along e 3 . In the transverse plane, the 
tensor T’ is locally 2-D and degenerate, therefore, it is identical to its isotropic tensor 
U’ and its deviator D’(also a local 2-D deviator of T) is singular. However, if the 
deviator D of T is singular then from Theorem 6, T has to be of triple degeneracy. 
Since T is only of double degeneracy, D must be nonsingular. 

On the other hand, if one of its local 2-D deviators D’ is singular, then the local 
2-D tensor T’ = U’ and T is degenerate. Since the deviator D is nonsingular, from 
Theorem 6, T is not of triple degeneracy. Therefore, T must be of double degeneracy. 

This completes the proof. O 

Singular points are the basic constituents of both vector fields and deviators. A 
singular deviator can be related to a vector with three zero components. It is natural 
to see the resemblance between them in topology. Since a tensor and it’s deviator have 
the same set of eigenvectors and their degenerate points occur at the same place, they 
share the same topological structure. This is why the topological structure of a first- 
order tensor(vector) and a second-order tensor field have such a striking similarity. 
In fact, although a second-order field is more complex than a first-order field, because 
the eigenvector fields of a second-order tensor field has sign indeterminacy, it usually 
has an even simpler structure. 


4.7 Representation of singularities in de viators 

The existence and location of degenerate points in a tensor field is determined by the 
Control Function. Analysis of the Control Function behavior in the neighborhood of 
a degenerate point also contains information regarding the distribution of singular 
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points in a tensor field especially in the case that singular points appear along a 
continuous line or surface. 

The control function A ( x , y, z) has maxima at a degenerate point xq, 


A (x 0 ) = 0 

dA ( j o) - n 
U 


(4.20) 


In the vicinity of this point, the control function can be expressed as: 

A(i) = A + (*-*„) 

t=l 9X i 

if 3 g \ 2 

( x i - x io) g^J A (so) + • • • (4.21) 

If a singular point is located on a continuous line or surface of singular points then 
the condition: A (x) = 0 is satisfied along the line or on the surface, respectively. 
Substituting Equation 4.20 into Equation 4.21, results in the following: 

A ( x ) = ^ ( x i - x m) A {x 0 ) 

+ I (| (ii ~ x<o) £:) AW)+ - 

= 0 (4.22) 


In Equation 4.22, the control function A (x) is expanded to the nth order so that 
at least one of its nth partial derivatives at x 0 is nonzero [76] [77]. Rearranging the 
function according to the order of variable (x — Xo)(here x — xq is written as x — x 0 , 
y — Vo and z — z 0 ), results in an implicit equation in 3-D space: 


o — fn,o,o(.x 3-o) n T fn—l,lfi(.X Xq) h 1 (y — y 0 ) + 


/n— 1,0,1 (x ~ X 0 ) n X {z — Z 0 ) H + 

/o,n,o(?/ ~ yoT + fo,0,n{ z ~ z o)" 


(4.23) 
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Depending on the values of the coefficients s, this equation is a representation 
of an object in 3-D space; surfaces, lines and isolated points are all possibilities [76]. 
The following is an example of an analytical tensor: 


T = 


x + y-\ 


x-y 


x-y -{x + y) + \ 0 


V 


(4.24) 




0 0 z 

By computing the trace and decomposing the tensor into its deviator and the 
isotropic tensor: 

5 2 rv \ 


D = 


x+y- f-f x-y 


V 


and 


x-y 

0 


^ 3 ( z + !) 


■(x + y) + l - | 
0 


0 


0 

0 


1 , 2z 

3 + T / 


(4.25) 


\ 


0 

u= o J(z + 1) 0 (4-26) 

0 0 i(z + l) J 

One can find the coefficients Q and 72(Equation 4.6 and Equation 4.14) of the 
characteristic equation, 


V 


Q — 2 (^11 "b ^22 ~b D33 "b 2 D l2 + 2 D l3 + 2 D 23 ) 


(4.27) 


= 


Dll D12 Di 3 

Di2 D22 D23 

D 13 D 23 D 33 

2 


-(-I) 4H)' 


1 

z ~2 


(4.28) 
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And then find the degenerate points as follows: 

• Triple Degeneracy (triple roots): Q = R = 0 
Result: an isolated point: ( x,y,z ) = (5,5,5)- 

• Double Degeneracy(double roots): (y) = — (^) ^ 0 
Results: 

• a line: x = 5, y = \ 

• a surface: (2 — ^) + (x — |) = 2 ( z ~ 2) 

The line and surface described above actually include both cases of double degen- 
eracy: 

• A x = A 2 > A 3 : x = y = | and z < \ for the line; \J{x — \) 2 -I- (x — |) 2 = 
^{z — |) for the surface. 

• Ai > A2 = A 3 : x = y = \ and z > \ for the line; \J{x — |) 2 + (x - |) 2 = 
— ^(z — |) for the surface. 

Figure 4.7 shows the resulting degeneracy: the small red ball is the location of a 
triple degenerate point; the line and the surface in the top figure is for X\ = A2 > A3 
case and the line and the surface in the bottom figure is for Ai > A 2 = A 3 . 

4.8 Topological Structure of a Singular Point in a 
3D deviator 

4.8.1 The Separating Surface of a Singular Point 

For second order tensor fields, in most cases, the eigenvector fields in the vicinity 
of a degenerate point can be described in terms of three types of angular sectors: 
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in the vicinity of the points of triple degeneracy. They also indicate that a locus of 
points of double degeneracy (A 2 = A 3 ) connects the points of triple degeneracy. This 
is evident from the two trisector points that lie in the symmetry planes just below the 
points of triple degeneracy. The existence of the line of double degeneracy is further 
verified by noting the two points of double degeneracy in the skeleton of the medium 
hyperstreamlines (Figure 4.13). 


4.9 Applications 

The following section contains examples of application using the methods discussed in 
this chapter applied to problems of scientific interest. We use both texture mapping 
and hyperstreamlines for display 2 . 


Deformation Tensor 


In the case of incompressible flow the deformation tensor, defined as 

De f=pl + pL 

OXj OXi 

has a zero isotropic part and therefore is equal to its deviator. 

For compressible flow, the deformation tensor is composed of a deviator superim- 
posed on a non-zero isotropic tensor which represents the rate of expansion. Therefore, 
a deviator describes the topological structure for both incompressible and compress- 
ible flows. 

For a JotationaL flow, inside. the. vortex core, flow is pure rotational. Assume that 
the flow advances in the z-direction and rotates around the z-axis while the velocity 
within the vortex core area is 

(-uy,ujx, 0 ) 


2 For dipole techniques, see Chapter 1 
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where u> is the angular velocity. By definition, the deformation tensor Def(r < R) 
becomes singular; outside the vortex core, the velocity is 


0 ) 


and its deformation tensor Def(r > R) is: 


-x 2 + y 2 0 N 
— xy 0 
0 0 } 

Here r is the distance from a point to the center of the vortex core and R is the radius 
of the vortex core. It is virtually a 2-D tensor with major and minor eigenvalues 
having equal magnitude but opposite sign and the medium eigenvalue remains zero. 
The deformation tensor is discontinuous at r = R. The angles of separatrices are 
calculated by using the tensor in the neighborhood of the vortex core Def(r > R) [38], 
and it turns out that there is no real solution for the angles. This indicates that the 
major and minor eigenvector fields are a pair of loci in the transverse plane while the 
medium eigenvector follows the direction of the vortex core. Studies of the alignment 
between vorticity and eigenvectors of the strain-rate (deformation) tensor in numerical 
solutions of Navier-Stokes turbulence have shown that the two principal strains with 
the largest absolute values (major and minor eigenvectors) lie in the equatorial plane, 
and the vorticity is automatically aligned to the intermediate eigenvector. 3 

Figure 4.14 (top) and Figure 4.14 (bottom) show the texture of a flow past a 
wingtip for a major eigenvector field and a minor eigenvector field, respectively. These 
two eigenvectors . remain in. the.. transverse, plane- perpendicular to the vortex core. 
Images are taken from a slice of the transverse plane along the vortex core, and the 
two eigenvector fields form two loci as we might expect. Color encodes the magnitude 
of their associated eigenvalues. 


2T 


xy 

- x 2 + y 2 


3 for detailed information, please see [27] [26]. 
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suggested that the model can be improved by taking into account the variation in 
swirling speed. 


4.10 Chapter Summary 

Tensor analysis is a very challenging task due to its complexity. Simplification is a 
natual choice. This chapter introduces a decomposition procedure to break a tensor 
field into a deviator and an isotropic tensor. A deviator carries the essential informa- 
tion about the tensor field, and describes the same topological structure. An isotropic 
tensor provides extra information like pressure in a stress field or expansion of a com- 
pressible flow in a deformation tensor, it is uniform at every location in the field yet 
can be dominant compared to the deviator. In order to understand the nature of a 
tensor field, we primarily study its deviator. 

Tensor invariants are very important quantities in tensor analysis. The character- 
istic equation of a general tensor has three invariants P, Q and R while a deviator 
only has two Q' and R'. These invariants determine that hyperstreamlines only ex- 
ist in between two surfaces in P — Q — R space and two curves in Q' — R' space, 
respectively. 

The degeneracy of a deviator is also its singularity, which explains the similarity 
between vector and second-order tensor fields. Because of the extreme importance of 
degenerate points to tensor analysis, the conditions for the existence of degenerate 
points are presented here. A control function is derived from these conditions which 
enables us to predict. the shape of this point set in a 3-D tensor field. The study of 
a control function sheds light on the representation of singularities in a deviator as 
well as in the tensor field itself. 

The singularities in a tensor field often have a very important physical meaning. 
For a deformation tensor, singularities only occur either at a vortex core or in a 
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constant flow; for a stress tensor in solids, it is a stress free point. After being 
decomposed into a deviator and an isotropic tensor, a stress tensor in a viscous flow 
with high pressure reveals the anisotropy of the flow. We applied our methods to 
several physical problems and the results are very encouraging. 


Chapter 5 

Conclusions and Future Research 


Vector and tensor fields are multidimensional multivariate data sets and are very 
difficult to comprehend. The objective of visualization is to simplify the analysis yet 
still capture the key tensor properties. This dissertation discusses the advantages and 
methods for visualizing and analyzing vector and tensor fields using topology. 

Because of its close connections with differential equations, combinatorial topol- 
ogy is the desired tool for studying vector and tensor fields. It deals with geometry, 
thus is able to find the intrinsic properties — critical and degenerate points and their 
corresponding separatrices that define the frame work of the fields (topological skele- 
tons). In this dissertation, research efforts on visualization developments for vector 
and tensor fields are based on topology analysis. 


5.1 Contributions 

5.1.1 Feature Comparisons for Vector Fields 

Much research has been done on vector field visualization. However, virtually none is 
on quantitative comparisons. In this dissertation, we define ” energy” as a quantity to 
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describe the feature distribution of a vector field. Two parameters a and 0 determine 
the local topology of a critical point. In the new a-0 space, all the basic patterns of 
simple critical points are distributed on a unit circle. Feature distribution is therefore 
defined as a set of critical points together with their a and 0 parameters. 

In Clifford algebra, vector fields are represented as a multiplication of linear fields, 
each of which is characterized by their a and 0 parameters. Therefore, a complex 
Clifford vector field is equivalent to a set of linear vector fields. This representation 
greatly simplifies the analysis. 

The concept ’’work” describes the energy difference between vector fields or the 
amount of work required when transforming vector fields. The Earth Mover’s Distance 
(EMD) is a linear optimization problem, it computes the minimal amount of work 
needed to transform one distribution into the other by moving ’’distribution mass”. 
EMD here is used to compute the least amount of work necessary to convert one vector 
field into another. It is a quantitatively measure of similarities and dissimilarities 
between vector fields. 

Multidimensional Scaling (MDS) computes a configuration of points in a low- 
dimension Euclidean distance to match the distance in the higher dimension as close 
by as possible. A vector field with a set of pairs of a and 0 values corresponds to 
a set of points in a-0 space. MDS in this dissertation is used to sort through a 
list of vector fields and produce a two-dimensional map that shows the relationship 
(distance) between pairs of vector fields. 

Feature comparisons using EMD is very useful when navigating through a large 
data base looking for vector Jields. with .similar topology or comparing results from 
experiments and computer simulations. It is ideal for data compression of a large flow 
field, since only the number and types of critical points along with their corresponding 
a and 0 parameters are necessary to reconstruct the whole field. It can also be used 
to better quantify the changes in time varying data sets. 
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5.1.2 Decomposition of a Tensor Field 

A second-order tensor usually has nine independent variables, even a symmetrical one 
has six. Therefore, a set of data from a three-dimensional second-order tensor field 
contains an enormous amount of information. Simplification of data before analysis 
becomes necessary. 

In this dissertation, a decomposition of tensors into a deviator and an isotropic 
tensor is introduced and serves as a preprocessor. The eigen system of a deviator 
is parallel to that of its associated tensor. It carries the essential information of a 
tensor field. An isotropic tensor serves as a uniform bias. Tensor invariants of a 
deviator describe the existence conditions of hyperstreamlines. Degenerate points are 
basic constituents of a tensor field. A degenerate point is also a singular point in a 
deviator. A control function is obtained from the characteristic equation of a deviator 
and it determines the occurrences of singular points in a deviator. The representations 
of singular points, both double and triple, can be derived from the control function. 

These singularities can further be linked to important physical properties of the 
underlying physical phenomena. For example, we show that for a deformation tensor 
in a stationary flow', the singularities of its deviator actually represent the area of the 
vortex core in the field; for a stress tensor in solids, the singularities represent the 
area with no stress; for a Newtonian flow, compressible flow and incompressible flow 
as well as stress and deformation tensors share similar topological features due to the 
similarity of their deviators; for a viscous flow, removing the large, isotropic pressure 
contribution enhances dramatically the anisotropy due to viscosity. 

5.1.3 Topological Structure of 3D Tensor Fields 

Topological representations of tensor fields reveal the geometric pictures that repre- 
sent tensor data locally as w^ell as globally in a very' simple manner. This dissertation 
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studies the topological structure of 3D tensor fields, an issue which hasn’t been dealt 
with before. The classification of degenerate points in 3D tensor fields is extended 
from that of 2D. Locally, separating surfaces divide the field around 3D degenerate 
points into several building blocks which are the fundamental elements in tensor fields. 
The separating surfaces have a general structure as they could appear at various an- 
gles as compared to 2D tensors which are confined on a plane. Each of the surfaces 
are characterized by patterns similar to those of hyperbolic or parabolic sectors and is 
bounded by hyperstreamlines that are emanating from or terminated at the degener- 
ate point. Consequently, a point of triple degeneracy can be classified by the number 
and type of separating surfaces surrounding it. The trajectories on the surfaces are 
locally 2-D, while off the surfaces they are fully 3-D. The set of degenerate points and 
their connecting separatrices form the topological skeletons that depict the geometric 
structure of the field without visual clutter. 

5.1.4 Representations of Degenerate Points 

Previous research efforts on tensor field topology only focused on isolated degenerate 
points. However, in 3D tensor fields, points of double degeneracy often appear as 
continuous lines and surfaces. This dissertation designs a method that derives the 
fully 3D representation from a control function of a deviator which works for general 
cases including both triple and double degeneracy. 

5.2 Analysis Framework 

Figure 5.1 and Figure 5.2 are two flow charts that summarize analysis processes of 
the vector field comparisons and 3D tensor field topology, respectively. 

Given n vector fields in a data base for comparison, Figure 5.1 starts with two 
input vector fields vl and v2 and the Clifford analysis procedure computes the feature 


Conclusions and Future Research 


147 



Figure 5.1: A frame work for feature comparisons of vector fields. 



Figure 5.2: A frame work for tensor field analysis. 
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distributions of these two vector fields, respectively. After obtaining feature distri- 
butions fl and f2, the next step is to compare whether these two distributions have 
the same number of critical points. If N 1 > N 2, Nl — N 2 regular points 1 with 
a = 0 and (3 = 0 are added to f2; if Nl < N2, N2 — Nl regular point 1 with a = 0 
and (3 = 0 are added to fl. EMD is computed when the number of critical points in 
fl and f2 are even. Similarly, the EMDs between all the pairs of n vector fields are 
computed. Then, MDS is performed and the n vector fields are displayed on a 2D 
map. 

Figure 5.2 describes the procedure for tensor field analysis. First, a tensor field is 
given as an input for the analysis. The decomposition step is a preprosessor to remove 
the uniform bias and obtain the deviator that describes the essential information of 
the input tensor. Eigen system is computed and degenerate points are located in step 
3. Then, the degenerate points are analyzed and their local topological structure are 
extracted. The integration step generates hyperstreamlines that connect the set of 
degenerate points and thus forms a global topological skeleton. In the end, a simple 
and compact representation of this input tensor is displayed. 


5.3 Future Research 

Vector and tensor fields are very rich in information. Current research covers only a 
fraction of them. There are still a lot of interesting issues for further exploration. 

• Feature comparison for 3D vector fields 

Currently, the method for feature comparisons only deals with 2D vector fields. 
A 3D critical point behaves similarly as a 2D critical point on the planes spanned 
by the pairs of eigenvectors of its Jacobian matrix. Therefore, it is possible to 
apply the feature comparison method to the 3D vector fields in a similar manner. 


Extension of feature comparisons for tensor fields 

A tensor field is equivalent to its eigen system. Eigenvector fields are essentially 
bi-directional vector fields, therefore, it is feasible to extend feature comparison 
methods to tensor fields. 

Addition of other characteristics of vector fields as criteria for feature 
comparisons 

Feature distributions used for comparison are not restricted to only topological 
characteristics. Other properties of vector fields can be used as criteria as well 
depending on the applications. For example, the locations of critical points in 
a flow are sometimes also important and are required to be part of the feature 
distribution. How to put totally different features together as one feature distri- 
bution set and assign different weights to them in order to rank the importance 
of these features is a particularly challenging problem. 

Numerical solutions for locating double degeneracy 

Representations of double degeneracy derived from control functions are still 
analytical solutions. Efficient computational methods to locate continuous lines 
and surfaces of double degeneracy will become necessary when dealing with real 
physical problems, both simulated and experimental. 

Further classification for triple degenerate points 

Patterns of 3D degenerate points presented here probably do not exhaust all 
the possibilities, and further classification maybe necessary for the prediction 
of 3D topological structures. 

Application of tensor field analysis to real physical problems Degener- 
ate points have very important physical meanings. For example, in the Boussi- 
nesq problem, the triple degenerate points are stress free points and in rotational 
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flow, the vortex core area is also an area of triple degenerate points of the defor- 
mation tensor. It will be very interesting and useful to apply these important 
properties to the actual physical systems for sovling problems such as vortex 
core detection. 


Appendix A 

Related Publications 


The work presented in this dissertation has appeared or will appear in several publi- 
cations; this appendix contains a brief overview of the relevant articles. 

• Yingmei Lavin, Rajesh Batra and Lambertus Hesselink, “Feature comparisons 
of vector fields using Earth Mover’s Distance”, Will appear in Proc. IEEE Vi- 
sualization ’98 

• Yingmei Lavin, Rajesh Batra and Lambertus Hesselink, “Vector Field Com- 
parisons using Earth Mover’s Distance”, Will appear in Proc. SIGGRAPH’98, 
Sketches 

• Yingmei Lavin, Yuval Levy and Lambertus Hesselink, “Singularities in Nonuni- 
form Tensor Fields,” in Proc. IEEE Visualization ’97, pp. 59-66, CS Press, Los 
Alamitos, CA., October 1997. 

• Lambertus Hesselink, Yuval Levy and Yingmei Lavin, “The Topology of Sym- 
metric, Second-Order 3-D Tensor Fields,” IEEE Computer Graphics and Ap- 
plications, March 1997. Special issue on scientific visualization. 
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Related Publications 
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Yingmei Lavin, Yuval Levy and Lambertus Hesselink, “The Topology of Sym- 
metric Tensor Fields,” in Proc. IEEE Visualization ’96, Late Breaking Hot Top- 
ics, pp. 43-46, CS Press, Los Alamitos, CA., 1996. 

Yingmei Lavin, Yuval Levy and Lambertus Hesselink, “The Topology of Three- 
Dimensional Symmetric Tensor Fields,” in 13th AIAA Computational Fluid 
Dynamics Conference, 1997. 
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1 Introduction 


Prim View provides a graphical user interface (GUI) and a set of development 
libraries to be used for scientific visualization. Reminiscent to development 
packages such as Open Inventor (and the corresponding visualizer iwiew), 
Prim View provides a package for the scientific visualization community. Cur- 
rently, Prim View is available only on the Silicon Graphics Workstation running 
IRIX6.2 and above. The libraries are written in C++ and use Open GL, and 
ViewKit. 

This document contains information required to create new and existing 
primitives from a developer’s stand point. The data created from these libraries 
can then be outputted to a file and viewed using a GUI package (see PrimView 
User’s Guide). A knowledge of C++ is required when working with these li- 
braries. 


2 Primitive Class 

The Primitive class is the parent from which all primitives are created. It stores 
information that is useful to its children. Table 1 lists the data accessible to its 
children. 


ID 

Unique identification for Primitive type 

Color 

Byte array storing R,G,B color values range (0-255) 

translate 

x,y,z offset from primitive’s origin. Initially identity. 

rotation 

Angle, and vector. Initially angle = 0. 

dirty flag 

Change of attribute -+ reconstruct primitive. 

Name 

String specifying label of particular primitive. 


Table 1: Primitive Class Storage 


In addition there are several methods that a developer of new Primitive 
children should override. These will be discussed in detail in Section 6. First, we 
shall discuss how to incorporate the current primitives supplied in the primitive 
library into your application. 


3 Using the Primitive Library 

Supplied with this manual is a SGI compiled library libPrimitives.so that you 
will link with your applications. Also included are several example programs 
using each of the various primitives. It is best to use these examples in con- 
junction with the header files included with the primitive library as learning 
sources. 

There are 3 basics steps to using any Primitive. These steps are shown in 
Figure 1. 
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Figure 1: Steps for using the primitive class 

The first step is to instantiate the primitive. Almost all primitive construc- 
tors contain optional parameters which can be filled either during creation or 
at a later time. A common constructor parameter is color. If no parameters are 
specified, default values are selected for the various constructors. In the case 
of colors, we attempt to select different colors for each primitive. In Section 5, 
default primitive colors are listed. The second step is to actually store useful in- 
formation for the primitive to display. Many basic primitives just require simple 
information such as a position or radius. Others, such as surfaces may require 
several hundred triangle positions (and normals). The third step is to output 
the primitive to disk. All primitives support the insertion operator, «. And can 
be outputted to a single file. Primitives will be concatenated automatically. In 
fact, two individual primitive files can be appended together with no problems. 
Every primitive is autonomous and contains information to reconstruct itself. 
Once the primitives are saved to disk, the data can be viewed in Prim View. 

4 Example 

To demonstrate the three basics steps to creating a primitive, we provide a 
simple example of creating some text to be displayed in Prim View. Granted, 
one could create the text directly in Prim View; however, using the library' has 
the advantage of creating multiple lines of text all within the same primitive. 
This provides for more efficient viewing than if each line of text were to be 
created as a separate primitive. 

The first line Text t instantiates the primitive. t.addO then proceeds to 
step two, and stuffs the data. Finally, file « t writes the data to a file. The 
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#include <f stream. h> 

#include "Primitives/Text.h" 

main (void) 

Text t; 

t.purgeO; //Wipe out any default text. 

t. add ("Creating Text" ,0. 1, .7) ; 

t.add("is real" ,0.4, .50) ; 

t. add ("Easy ! " ,0.6, .30); 

of stream f ile("text .pv" ,ios: rout) ; 

file << t; 

f ile.closeO ; 

> 


Figure 2: Example creating a Text primitive 



Figure 3: Output from Text primitive code (see Figure 2) 
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output of Figure 2, which can be viewed directly from Prim View is shown in 
Figure 3. Other examples are provided with the library, and follow the basic 
three steps. 


5 Primitives provided 

Currently, there are 12 basic primitives provided within this library release. 
Below is a description of each primitive and useful notes to keep in mind as you 
use them. For detailed use of the primitive see the respective header (.h) file. 

Primitive: Base class which all other primitives must inherit from. Discussed 
in detail in Section 6. 

Axes: Cartesian coordinate system with x,y,z labels. Default color is white 
(R=G=B=255). The coordinate system can be moved around via the 
origin () attribute, and the coordinate system can be scaled using the 
size() attribute. Both attributes can be called without parameters for 
querying. 

Arrow: Creates an arrow with full control over its characteristics. Characteris- 
tics include head (head is where the arrow head is located) and tail location 
(head, tail). The size of the head radius which is the largest part of the 
head (it tapers to a point) can be controlled via headRadius(). The tail 
radius, which is the radius of the barrel is controlled using tailRadiusQ. 

The location of where the head of the arrow starts relative to the en- 
tire arrow is specified with (headPercentageO using a value between 0.0 
to 1.0 corresponding to a percentage. The arrow’s default color is green 
(0,255,0). 

Cube: A Cube is composed of 8 points. They are specified individually us- 
ing the addElemO method. See header file for description of order in 
which cube is connected. The default color for the cube is red (255,0,0). 

The line thickness (glLineWidth in Open GL) default is 5.0. The point 
size (glPointSize) is 2.0. These can be changed via setLineSizeO , 
setPointSize () . The point color is set using the standard method setColor () . 
The line color is set using a new method setLineColor(). 

Grid: The Grid class is a wrapper for a plot3D structured grid. Providing the 
grid class with a Plot3Dio class (p3dio()) will cause the Grid to import 
in a structured grid. The color by default is olive green (0,127,0), and will 
be drawn as a mesh with k=0 shell selected (i and j is drawn to the full 
span). The step size for each i j,k slice is 1. The line thickness is default to 
1.0. The grid can be rendered solid using drawStyle() . The line thickness 
used to draw the mesh is controlled via lineSizeO . There are times when 
only a subsection of the grid needs to be drawn or calculated. iMin/ iMax , 
jMin/jMax, kMin/kMax are attributes that control the scope of the grid. 
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IsoSurf: Designed to work in conjunction with algorithms such as marching 
cubes. Triangles are drawn in the order specified using addElem(). The 
default color for the surface is gray (95,95,95). In the constructor, a hint 
can be given to the IsoSurf object as to how many triangles are expected 
to be inputted. Internally, memory is allocated to the hint size. When 
memory has been exhausted by specifying more triangles than the hint, 
memory is reallocated to double the hint. This continues until all triangles 
can fit into memory. Reallocation is an expensive operation, therefore a 
proper hint can greatly increase efficiency. 

Legend: Given a ColorMap is all that is needed to construct a legend. The 
ColorMap controls color space (RGB/HSV), and the minimum and maxi- 
mum values to interpolate between. The Legend can be positioned using 
the setOriginO method. The origin corresponds to the bottom left hand 
comer of the Legend. The coordinates are specified in 2-D and correspond 
to a percentage of the screen (scaled between 0.0 to 1.0). The coordinates 
follow a first quadrant system (i.e. origin at bottom left hand portion of 
the screen.) 

Line: Lines are composed of a series of points. The points are added in linear 
fashion using the addElemC) method. The points can be changed later 
using the setOriginO method. (To get total number of points added 
to Line use the Length () attribute.) Lines can be drawn as a series of 
points and not connected (use set Connect ion () . Points can have a differ- 
ent color than the line itself and different point sizes (setPointSizeO). 
The color of the line and its line thickness can be controlled as well 
(setLineColorO , setLineSizeO). As a convenience function, both the 
line and point color can be changed simultaneously (setColorO). Both 
point-size and line-size are set to 3.0 by default. Line color is an off blue 
(purplish) (165,0,255). 

Node: Similar to Line, however, Node can specify only one point. Nodes 
can be drawn as a simple Open GL point or can be rendered as a sphere 
(representation ()). Default point size is set to 3.0. Default sphere 
radius is set to a tenth of the point size. 

TexSurf: is used to render textured surfaces. The items needed for a texture 
are a Texture (setTextureO), a corresponding surface to map texture 
to (setSurf ace()) and a scalar value map to set each textured pixel 
(setColorValues () , setColorMapO). TexSurf should be used in con- 
junction with the TextureCalc class which is a texture calculator designed 
to create arbitrary surfaces and compute textures (via vector fields) over 
these surfaces. 

Text: Designed to create text annotations. This class can add an arbitrary 
number of text elements at arbitrary locations via the add() method. 
Coordinate locations are specified in 2-D and are a percentage of the 
screen coordinates (between 0.0 to 1.0). The bottom left hand corner is 


6 


the origin. Text sentences are referenced from the bottom left hand corner 
as well. 

Tube: The tube class is a general class creating arbitrary tubes in 3-D space. 
A tube is created from a series of ellipses each possessing a color. The 
tube merely stitches the ellipses together and linearly interpolates the 
color between ellipses. Ellipses are added in sequential fashion using the 
addElemO method. The resolution of the ellipses can be controlled af- 
ter the ellipses have already been added (setEllipsesResolutionO). 
Ellipses contain a scalar value which is interpolated and mapped to a 
color using ColorMapO. The color map can be changed at any time. 
In addition, the ellipses’ scalar values can be overridden by setting the 
colorStyleO attribute to UNIFORM. The color specified via setColorO 
will be used. 


6 Creating Primitives 

All new primitives created by the user must inherit from the base Primitive 
class. This parent class provides header information (on a primitive by primitive 
bases) that makes your primitive recognizable throughout the system. 

6.1 Constructor 

When creating a new class, a unique primitive ID needs to be registered for the 
primitive. In the file PrimitivelD.h is a table of enumerations. New primitives 
must be appended to this table. It is recommended that new primitives from a 

particular group/agency start at a unique base i.e. enum ePrimitivelD { . . . TEXT, DOD_SURFACE=200}. 
If the primitive is useful and available for others, it can be incorporated into 
the base library and distributed as such. 

The remaining parameters in the Primitive constructor are as follows: 

Primitive(ePrimitiveID=PRIMITIVE, const char *name=NULL, 

unsigned char R=95, unsigned char G=95, unsigned char B=95) ; 

The parameters all have non-essential, but highly useful data elements. A 
name is useful for identifying the primitive to a user. This name is what is 
shown in PrimView (member data _primName). The color parameters should be 
used to specify the color of the primitive (member data .color [3]). 

6.2 I/O and versioning 

When users invoke the insertion and extraction operators (»,«), calls to 
readlnO, vriteOut() are made. By default, the Primitive’s readlnO and 
writeOutO are invoked. The data written/read by this default are the primi- 
tive ID, the version number of this primitive class, the name (in binary), color, 
and local rotation and translation vectors. The version number is acquired by 
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void Primitive: :setColor (unsigned 

unsigned 

unsigned 

{ 

_ color [0] = R; 

_color[l] = G; 

_color[2] = B; 

_dirty = 1; 
notify () ; 


> 


char 

char 

char 


R, 

G, 

B) 


Figure 4: Example of using notify and the -dirty flag 


invoking the method version (). Hence, developers are encouraged to override 
this attribute to return a version for their primitive. In addition, another at- 
tribute oldestVersionO should be overridden to support the oldest version 
that can be loaded by the class. Ideally, we would like backward compatibility. 
Using the oldestVersionO method is used by the default readlnO method 
to determine whether the primitive can be loaded into the system. If it can’t, 
the read process is aborted with istream’s fail bit set. If the version information 
is valid, the readlnfistream &in, int version) method is invoked. The de- 
veloper should then use the version information to load the appropriate values 
(and set default values for information not specified in an older version). 

6.3 Use of the dirty flag and Notification 

The primitive class uses a simple notification system. Other classes using the 
primitives can register itself to be notified any time a change is made within 
a primitive. Therefore, in order to correctly reflect any changes made within 
a primitive (typically screen redraws) it is important to invoke the notifyO 
method. To clarify this point let’s look at the setColorO method shown in 
Figure 4. 

Once the color has been stored, the -dirty flag is set to true, and the function 
notifyO is invoked. The notifyO method informs any registered users that a 
change has been made to the primitive. A typical registered user is PrimView’s 
View engine. Once the View is notified that a change has been made to the 
primitive, it invokes the primitive’s drawThySelf () method. The -dirty flag is 
internally used for efficiency. For example, in the drawThySelf () method, the 
typical Open GL construct is to use call fists. This greatly increases efficiency 
by invoking the call fist instead of recomputing the primitive. However, when a 
change has been made to the primitive (color, radius, etc.) it maybe necessary 
to recompute the primitive. Using the -dirty flag informs the drawThySelf () 
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method that recomputation is necessary. Of course, once the computation is 
complete, the _dirty flag should be set to false. 

6.4 Drawing the primitive in Open GL 

Since the drawThySelf () method uses Open GL to draw itself - developers 
should keep the template shown in Figure 5 in mind. 
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void YourClass: rdravThySelf (void) 

{ 

//Enter a local coordinate system 
glPushMatrixO ; 

//Apply any local coordinate transformation 
glRotatef (.rotate [0] ,_rotate[l] , _rotate[2] , .rotate [3]) ; 
glTranslatef (.translate [0] ..translated] , .translate [2] ) ; 

//Check if dirty flag set. If so recompute primitive, 
if (.dirty) 

{ 

//Clear the dirty flag 
.dirty = 0; 

//Use Open GL commands to create a call list 
//.iDisplay contains the old display list ID. 
if (.iDisplay) glDeleteLists (.iDisplay, 1) ; 
.iDisplay = glGenLists(l) ; 

glNewList (.iDisplay, GL.COMPILE.AND.EXECUTE) ; 

{ 

// YOUR CODE HERE. 

} 

glEndListO ; 

} 

else 

{ 

glCallList (.iDisplay) ; 

> 

glPopMatrixO ; 

> 


Figure 5: Template for drawThySelf () method 
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1 Introduction 


Prim View provides a graphical user interface (GUI) and a set of development li- 
braries catered for scientific visualization. Reminiscent to development packages 
such as Open Inventor (and the corresponding visualizer ivview ), PrimView pro- 
vides a package for the scientific visualization community. Currently, Prim View 
is available only on the Silicon Graphics Workstation and runs on IRIX6.2. The 
libraries are written in C++ and exploit Open GL, and ViewKit. 

This document contains information required to run the Graphical User 
Interface portion of PrimView. Users wishing to learn how to create data to be 
viewed by the GUI should refer to the library guide . 

2 Starting Prim, View 

PrimView can be started with or without a command line argument. Starting 
without an argument loads the basic GUI interface with no data. An optional 
argument can be given which specifies a PrimView file (typically denoted with 
a .pv file extension). 

On startup, the camera (which is what the user is looking through) points 
down the negative z-axis. The camera is situated a fixed percentage away from 
the bounding sphere of the loaded primitive. As additional primitives are added 
to the world (or scene) the bounding sphere is recomputed. However, the camera 
is not reset if the world contains at least one primitive. To reset the camera to 
contain all elements in the world see the pan function in Section 3. As primitives 
are added a bounded box primitive is continually updated to encompass all 
primitives in the world. This primitive is normally hidden but can be displayed 
(see Section 5 for displaying the bounding box). 

3 Prim View Window Components 

Figure 1 shows the major areas of the PrimView window. A brief explanation 
of each region follows: 

View Provides a depiction of the primitives contained within the world. 

Left mouse button simulates a trackball and rotates the view about its 
own local coordinates. 

Middle mouse button pans or translates the world. 

Left and Middle button Holding both the left and the middle mouse 
button and moving the mouse up causes the world to zoom out. 
Sliding the mouse down while holding both buttons causes the world 
to zoom in. Note that each region is context sensitive and therefore, 
all mouse clicks must be within the view window for these actions to 
occur. 
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Figure 1: PrimView window 


Right Mouse Button Provides options that manipulate the view. These 
menu choices are elaborated in Section 4. 

Console Provides transformations on the View region. 

X Y Z spindles Rotates world (located in View) about screen X,Y, and 
Z axes. See View for local rotations. 

Pan buttons Translates world in screen X, and Y coordinates. Pressing 
the button labelled Pan causes the camera to reset the view. 

Zoom In/Out Zooms the world by manipulating the camera. 

Options Controls which primitives are shown in the View, and allows for 

editing of individual primitives. 

Primitive Namelist is the region that contains the names of all of the 
primitives loaded in the world. Names preceded by an asterisk () 
indicate that the primitive is hidden and will not appear in the view. 

Filter List Menu filters the list of items depicted in the Primitive 
Namelist region by primitive type. Selecting All shows all prim- 
itives in the Primitive Namelist region. 

Right Mouse button Manipulates the primitives within the world. See 
Section 5 for details. 
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4 View Region Menu Options 

Pressing the right mouse button over the View region yields a menu that pro- 
vides manipulation of the View and the world (scene). Following is a description 

of each menu item found in the View menu with a brief description. 

New... purges all primitives from the world, and resets the view as described 
in Section 2. 

Open... adds primitives found in primitive file to existing world. If the primi- 
tives are older versions, Prim View will upgrade them to the most current 
version. Camera view remains unchanged as primitives are added.. 

Import FAST grid... This command will detect and load four types of Struc- 
tured grids. Single-Zone, Single-Zone with I-Blank, Multi-zone (struc- 
tured), and Multi-zone (structured) with I-Blank. The imported grid is 
converted to a Grid primitive which allows additional attributes to be 
changed and saved (name, color etc.) 

Save Saves all visible primitives in the view region to disk. The filename used 
under the saveas... menu option is used. If no filename has been specified, 
the saveas... option is automatically invoked. 

Save As... Saves all visible primitives in the view region. A file requester 
prompts for the filename. 

Save RGB... Saves current view region as an SGI RGB picture file. The size 
of the RGB file is exactly that of the view region. 

Save View... The camera angle, position, frustum and background color set- 
tings are saved to disk. This allows particular camera views to applied to 
the view region. 

Load View... Loads the camera, frustum and background color settings and 
applies it to the view. 

Background Color Brings up a color wheel, allowing the background color of 
the view region to be changed. 

Reset View Has the same affect as pressing the Pan button in the console. 
Resets the camera such that the camera is placed along the z-axis and is 
pointing in the negative direction. 

Options Toggles the visibility of the Option Region (see Figure 1). 

Show Console Toggles the visibility of the Console Region (see Figure 1). 

Full View Toggles between hiding both the Console and Option Region to the 
current Region view. 
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Perspective When toggled on (yellow box visible) perspective view is on (fore- 
shortening of lines), else Orthographic projection is set. 

Lighting Lighting produces shading of primitives. Without lighting pixels col- 
ors are fully saturated and uniform. 

Two Sided Lighting Primitives whose normals are not pointing in the same 
direction as the light yield dark regions. Two Sided Lighting defeats this 
by lighting both sides equally. In effect as if to put a light 180 degrees to 
current light. 


5 Option Region Menu Options 

All primitives that are in the world (scene) can be controlled via the Option 
Menu. Individual primitives can be selected by holding down the control key 
and pressing the left mouse button over each primitive name located in the scroll 
window. Pressing the right mouse button drops a list of options for manipulating 
the selected items. 

5.1 Hide, Show, Remove... 

Hide hides primitives from the View Region. Primitives are still in memory 
and are denoted by an a prefixed asterisk (*) in their names shown in the 
Option scroll window. 

Show shows all selected hidden primitives in the View Region. 

Remove... On confirmation, deletes selected primitives from world (freeing up 
memory). 

5.2 Edit... 

PrimView provides some basic editing capabilities on the primitives loaded 
within the world (scene). Elements that have been selected can be edited by 
selecting the Edit... menu option. In many cases, selecting multiple primitives 
allows the changing of each primitive’s attributes collectively. For example, to 
change the name of all arrows to ’Force Vectors’, one can select the Arrow op- 
tion in the filter menu, select all arrows with the Control key, and then select 
Edit from the Option Menu. Primitives not of the same type can not be edited 
together. 

By in large, the majority of the primitives allow their name and color to be 
changed. Exceptions are noted below. In the Edit form, the Okay button accepts 
all changed and closes the Edit form. The Apply button merely previews the 
changes. To make the changes permanent the Okay button needs to be pressed. 
To cancel the changes made by the Apply button, press the Cancel button. This 
too causes the Edit form window to close. 
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Below is a list of the various primitives that can be edited and notes to help 

clarify the edit form. 

Arrow : Arrows are used typically to identify locations of interest. The Edit 
form allows the arrow’s head and tail to be manipulated similar to a vector. 
In addition, the Arrow Head Radius can be changed as well as the Tail’s 
radius. The starting point of the Head can be specified as a percentage 
(between 0.0 (no head) to 1.0 (all head) ) 

Axes The position of the Axis and the size can be controlled via this form. 
The size is useful for scaling purposes. 

Cube Not implemented in this version. 

Grid Allows dynamic control over which grid slice is currently displayed. This 
is controlled via the min/max sliders and the corresponding selected grid 
slice toggle. For example, to control the ’k’ slice of a grid, make sure the 
’k’ toggle is selected and slide the min or max slider. On low performance 
machines, the min and max can be inputted directly in the form and either 
the Apply or Okay button can be pressed for the change to take place. 
The Step controls the coarseness of the grid. For example, selecting a step 
size of 2 would draw every other point to form the grid. This is useful for 
low performance machines. A grid can be displayed either as a mesh or 
as a shaded surface. This is controlled by the View As toggles. 

Isosurface Currently the name and color are the only properties that can be 
modified. 

Legend The Legend is considered a separate primitive and is not linked to the 
primitive that spawned it. The position and scale size can be controlled. 
The X and Y value vary between 0 to 1 and correspond to screen position 
of the legend’s bottom left comer. (0,0) corresponds to the bottom left 
hand comer of the screen and (1,1) is the upper right hand comer. 

Line Currently not implemented. 

Node Nodes can be represented as points or as spheres. In either case the 
position and scale size can be modified. 

Texture Color cannot be modified. Depending on the capabilities of the dis- 
play machine hardware or software texture mapping can be selected. For 
systems supporting hardware texture mapping set the toggle to Texture 
Mapping. For systems that do not have hardware texture mapping, tex- 
tures are rendered using triangles when the Triangle Rendering toggle is 
selected. Note: even on platforms which do not support hardware texture 
mapping, OpenGL will internally render the texture in software if the Tex- 
ture Mapping toggle is set. Tests have indicated that Triangle Rendering 
has better performance than OpenGL’s software texture mapping. 
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Text The position of the text’s bottom left hand corner can be controlled in 
screen coordinates (0,0 - bottom left hand corner of screen 1,1 - top right 
hand comer of screen) via the Origin. The Content field changes the text 
displayed in the View region. The Name is merely a tag used to identify 
the text in the Option’s Scroll window. 

Tube Tubes can either be single colored or by representing a scalar value at each 
ellipse can depict multiple colors. Colors are formulated by interpolating 
within either an RGB or HSV colorspace. By default the minimum value 
is assigned blue and the maximum value assigned magenta. Colors in both 
spaces span the following spectrum: Blue, Cyan, Green, Yellow, Red, and 
Magenta. Under the Color Control field a tube’s maximum and minimum 
value can be set. By changing this range, the corresponding color space 
interpolation changes. Colors used in the spectrum can be modified by 
selecting the spectrum option located in the color control region. Note that 
the spectrum must maintain its order i.e. a tube can span Green (min) 
to Yellow (max) but never Yellow (min) to Green (max). Selecting HSV 
or RGB selects which color space to interpolate in. RGB typically yields 
a much more discrete space best for scientific visualization. HSV tends to 
smooth the colors into one another and typically yields a more aesthetic 
representation. Selecting Single Color... overrides the tubes scalar data 
and displays the tube as a single color. The scalar data is not lost, and 
can be recovered by pressing the Unify Tubes... button. Unify Tubes... is 
typically used when multiple tubes have been selected. Its function is to set 
all of the selected tubes to a single minimum and maximum corresponding 
to the selected tubes’ overall minimum and maximum values. This is useful 
when a uniform color map is needed over a set of tubes. Spawn Legend will 
create a legend primitive based on the minimum and maximum value as 
well as the color space found within the form. A tube’s corresponds to the 
number of points used to draw each ellipse. The scale affects the overall 
radius of the tube. Note: Anytime the scale is changed, the value reverts 
to 1. This indicates that the current scale is now 1 and any additional 
change will affect the current tube’s scale. 

5.3 Create > 

Certain basic primitives can be created directly within Prim View. These prim- 
itives are typically used to provide supplementary information to View Region. 
Currently, there are four primitives that can be created from the Option Menu: 
Arrow, Axes, Node, Text. Creation of any of these primitives automatically 
nvokes the edit primitive form. 
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1 Abstract 

A method for comparing three-dimensional vector fields 
constructed from simple critical points is described. This 
method is a natural extension of the previous work [1] which 
defined a distance metric for comparing two-dimensional 
fields. 

The extension to three-dimensions follows the path of 
our previous work, rethinking the representation of a critical 
point signature and the distance measure between the points. 
Since the method relies on topologically based information, 
problems such as grid matching and vector alignment which 
often complicate other comparison techniques are avoided. 
In addition, since only feature information is used to repre- 
sent, and therefore stored for each field, a significant amount 
of compression occurs. 

2 Introduction 

Vector fields 1 are used to study phenomena in almost all 
areas of the physical sciences including such diverse subjects 
as climate modeling, dynamical systems, electromagnetism, 
and fluid mechanics. Hence, the analysis of vector fields has 
become a significant concern to the sciences, and a variety of 
techniques for visualizing and analyzing vector fields have 
been developed. However, an effective technique to quanti- 
tatively compare vector fields has not been developed. This 
paper addresses the issue. 

A review of existing comparison techniques is first dis- 
cussed. Even the most promising of these techniques lack 
the quantitative capabilities for automated comparisons. The 
properties of the two dimensional classification technique 
used in [1] are briefly discussed and are used to extend the 
classification to three-dimensional critical points. A com- 
plete categorization of 3-D simple critical points is presented 
and is used to redefine the EMD metric allowing for a quan- 
titative comparison between 3-D flow fields. The paper con- 
cludes with an example demonstrating the effectiveness of 

' The definition of vector fields is restricted to continuous fields or flows 
which are discussed in section 4. 


the technique on a thermal convection model described by 
the Lorenz equations. 


3 Existing Comparison Techniques 

A variety of comparison techniques exist for vector fields. 
These techniques basically fall into three general categories: 
Image, data, and feature extraction based comparisons. In 
most of these cases, comparisons are made visually [2], 

Image based comparisons work on the computer gener- 
ated image. Often times, a numerical data set is converted 
into an image that simulates an experimental visualization 
technique (computational flow imaging). This may be eas- 
ier than extracting a vector field from an image, such as 
Schlieren. However, visualizing a field in 3-D is quite diffi- 
cult and often, these techniques are limited to two dimen- 
sions. In addition to side-by-side comparison of images, 
other techniques include image fusion, and Fourier analy- 
sis [3], 

Data level comparison techniques operate directly on the 
raw data. An accurate comparison requires proper grid 
alignment which can involve problematic interpolation be- 
tween two fields [4]. 

The last comparison category is the extraction of features. 
Typically features are problem specific; for example in fluid 
mechanics features include vortex cores and shock surfaces. 
Often there is a geometric representation of the feature and 
possibly a semantic representation of the system which can 
be compared using a pattern recognition technique [5]. This 
may lead to more robust comparisons. 

Qualitative comparisons have been based on the concept 
of critical points in vector fields. Past study has focused on 
the geometric structure of vector fields [6] and last year a 
quantitative measure for two dimensional vector fields was 
introduced [1]. The work is extended by defining a quan- 
titative measure of the similarities and differences of three- 
dimensional vector fields. 


4 Description of Phase Portraits with 2-D a-/3 
Parameters 

Since the method reduces a 3-D flow pattern into 2-D 
components, the relevant 2-D categorization for this type of 
vector field is reviewed. A 2-D vector field that can be rep- 
resented as a system of two simultaneous differential equa- 
tions has the following form: 

Vx = ^: = F(x,y) (1) 

Vy = ~j~ = ^ 

where F and G are continuous and have continuous partial 
derivatives in some region D. The solutions to this system 
forms a family of directed paths. Given some initial value 
to the system, a parametric representation expressed as x = 
<f>(t), y = ip(t) can be deduced. The image formed is the 
phase portrait and is typically described by the number, type, 
and arrangement of critical points (or equilibrium points). 
These are points where F(x,y) = 0 and G(x,y) = 0. The 
nature of a critical point will not change under continuous 
(affine) transformation. Critical points are significant in that 
they are the only points in a vector field where tangent curves 
may cross each other. Therefore, critical points delineate the 
field into sectors of uniform flow. 

A critical point is said to be isolated or simple if there is 
an open neighborhood around it that contains no other criti- 
cal points. The behavior of the flow about a critical point can 
be analyzed by investigating the trajectories in the neighbor- 
hood of the critical point. If the distance is sufficiently small 
(say dx, dy), a first order approximation (Equation 2) of the 
field can be used. 

dv dv 

v x (dx, dy) » ~gf dx + ~Q^ d y ( 2 ) 

v y (dx, dy) « ~^dx + 

Hence, the flow pattern is completely determined by the Ja- 
cobian, Jij = (i,j = 1,2) evaluated at the critical point. 
The various patterns formed in the phase-plane space can be 
seen by analyzing the eigenvalues of the Jacobian. The char- 
acteristic equation 


A 2 + FA + Q = 0 (3) 

where P = —trace(J) and Q — det(J) is used to classify 
the various patterns using the well known P — Q stability di- 
agram [7]. However, advantageous properties arise by defin- 
ing a new space {a',/3') as explained in [1], wher e the eigen- 
values map a = P and (3 — sign(P 2 — 4Q)y/\P 2 — 4<2| 
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Figure 1 : Basic patterns for simple critical points. 


and are normalized as follows: 


a 


& 


s/d 2 + 13 2 
(3 

s/a 2 + 3 1 


(4) 

(5) 


In this space, critical points obey all the the rules defined for 
a regular 2-D Euclidean space and the distance between any 
two critical points is a metric. 

It is shown in [8] that the actual values of a and /3 do 
not determine the portrait of the critical point; only the ra- 
tio between them matters. Hence, this normalization maps 
all points onto a unit circle and thereby provides a means of 
relatively quantifying the difference between various points 
by just an angle. Also note that a uniform vector field with 
no critical points, F(x,y) = ci, G(x,y) — ci has a — 0 
and j3 = 0 and maps to the origin of the unit circle. This 
is the reason why arc length is not used as a metric. For 
the remainder of the paper, a and j3 values will be assumed 
normalized. The patterns are sketched in Figure 1 and enu- 
merated in Table 1 2 . Notice a positive or negative real part 

2 The definition of saddle indicated in the table is more relaxed than spec- 
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Type 

Constraint 

> 0 

= 0 

Repelling Star 



> 0 

> 0 

Repelling Node 

101 < 

M 

— 

> 0 

Saddle 

\8\ > 

jaj 

< 0 

> 0 

Attracting Node 

101 < 

|a| 

< 0 

= 0 

Attracting Star 



< 0 

< 0 

Attracting Focus 



= 0 

< 0 

Center 



> 0 

< 0 

Repelling Focus 




Table 1 : Classification of Critical Points via a- 8 values 

(denoted by a) is indicative of repelling/attracting behav- 
ior. And if an eigenvalue has an imaginary part (/? < 0), 
it indicates circulation about the point and the trajectories 
can be represented via logarithmic spirals, otherwise asymp- 
totic behavior whose trajectories can be described via simple 
power laws is exhibited. 

5 Classification of Three-Dimensional 
Vector Fields Using Phase Portraits 

The formulation for a 3-D vector field is very similar to 
the 2-D analysis. For a three-dimensional vector field, the 
Jacobian is represented by a 3X3 matrix, Jy = (i,j = 
1, 2, 3) The characteristic equation now becomes 

A 3 + PA 2 + QX + R = 0 (6) 

where P = —trace(J), Q = \{P 2 — trace(J 2 )), and 
R = —det(J). Three distinct eigenvalues are possible, 
along with three eigenvectors. The flow field can be de- 
composed into fundamental solution trajectories along its 
eigenvector planes as demonstrated by Reyn [9] and Chong 
et al. [10]. All other solutions trajectories converge (or di- 
verge) to these eigenvector planes. Therefore a critical point 
in 3-D can be defined by a set of three (a, 8) values. Each 
(q,/ 3) point corresponds to a solution trajectory formed in 
the respective eigenvector plane. 

To simplify the classification of the various phase por- 
traits about a three-dimensional critical point, the Jacobian 
is transformed into canonical form. This does not affect 
the eigenvalues since they are invariant to changes in scale, 
translation, and rotation. Philippou, and Strickland catego- 
rized the Jacobian into seven basic canonical forms [11] or 
classes. With each form, several phase-portraits are possi- 
ble. In Tables 2 and 3, all possible cases are enumerated in a 
similar style as presented in reference [1 1] 3 along with the 

ified in the figure. 

3 The class structure is slightly changed by placing the complex Jordan 

form last. 



Figure 2: Decomposition of a class 1 critical point along 
eigenvector planes. 


corresponding a, 8 sets. A brief discussion of the various 
classes is given below. 

Class 1 (Ai, A 2 ,A 3 ): Class 1 is indicative of eigenval- 
ues which are real and distinct In this case, there exists 
three independent eigenvectors and therefore three indepen- 
dent eigenvector planes. For the case of a Hermitian matrix 
(and its subclass the real symmetric matrix), the eigenvectors 
are mutually orthogonal. If all of the eigenvalues are posi- 
tive, repelling nodes form in all three planes. If the signs 
differ, a saddle occurs in two of the planes, and finally if 
all signs are negative attracting nodes occur. All other solu- 
tion trajectories approach or diverge from the critical point 
as t — > oo and are not planar. In this case, there exists a full 
set of a, 8 values. A degenerate 2-D case is exhibited when 
one of the eigenvalues is zero. Only one plane will contain 
a simple solution trajectory, the other two planes will con- 
tain lines (since in 2-D every plane in the third dimension 
is identical). In this case, there exists only one distinctive 
a, 8 value. The remaining two are set to zero. In fact, all 
2-D cases degenerate to {{ai,/3i} , {0,0} , {0,0}} and the 
computed comparison values are identical to those in refer- 
ence [1]. Figure 2 is an example of trajectories formed for a 
node/saddle/saddle combination. 

Class 2 (Ai , Ai , A 2 ): This is a degenerate case where two 
eigenvalues are identical. The multiplicity is 2, however, 3 
independent eigenvectors can still be found. One plane will 
contain a star pattern, the planes normal to this will contain 
solution trajectories (nodes or saddles) and will have identi- 
cal a, 8 values. The 2-D case degenerates to the star pattern. 
The 1 -D (Aj = 0) case is ignored. 

Class 3 (Ai,Ai,A 2 ): In this case, only 2 independent 
eigenvectors and only two solution trajectories exist In one 
plane, a log star 4 trajectory is observed and in the other plane 

4 A log star is also referred to as an improper node see [12] 


a node or saddle. A log star in the a, 0 space is indistin- 
guishable from a star pattern, i.e. a = ± 1 , 0 = 0 since the 
star pattern is just a special case of the family of logarithmic 
stars formed. 

Class 4 (Ai, Ai, Ai): Case 4 exhibits triple degeneracy 
with three independent eigenvectors. Any plane passing 
through the critical point will exhibit a star trajectory. The 
set of a, 0 values are identical. 

Class 5 (Ai , Ai , Ai ): Two linearly independent eigenvec- 
tors exist for this triple degeneracy case. Here there are only 
two independent planar trajectories. One trajectory is a log 
star located in the coordinate plane spanning the eigenvector 
(in x-y plane for the canonical form) and a star pattern in the 
other coordinate plane (x-z). 

Class 6 (Ai, Ai,Ai): In this case, the multiplicity is three 
but only one independent eigenvector exists. Therefore, only 
one plane contains an attracting/repelling log star trajectory. 
Hence only one unique a, 0 value exists and this case con- 
flicts with the 2-D case of class 2 and 3. Therefore, some 
false positives can be expected. Fortunately, class 1 is the 
most common occurrence [ 11 ], which this method classifies 
most uniquely. 

Class 7 (a + 0, a — 0, A 3 ): Three eigenvalues are 

found, two of which must be complex conjugates of each 
other for the J matrix which contains no imaginary values. 
Only one plane will contain solution trajectories which are 
either a focus or a center. Hence, only one unique a, 0 pair 
value exists just as in the 2-D case. The real eigenvalue, 
A 3 , denotes a stretching or compressing phenomenon where 
trajectories either spiral away or towards the solution plane, 
and this will not be captured [13]. 


associated with the vector field’s critical points: 


{{{a?\0[\(a£\0<\(4\0P)}, 

{(a?\0i\(a?\0i 2 \(a?\0™)}, 

... ,{(a[ n \0[ n) )0a^\0 { 2 n] )0ai n \0i n) )}} ( 7 ) 

Definition 2 (Energy) The energy for a vector field is: 


Energy = , I ££<(«“>)’ + <#>’>, 
j=l i= 1 

where n is the total number of critical points in this field. 


The energy is a quantity that characterizes the critical points 
of a vector field. It is different from physical energy. The 
concept “work” is used to measure the energy differences 
between two vector fields or the amount of energy used to 
transform one vector field into the other. For a 3-D vector 
field, work will be defined at two levels. At the higher level, 
the work required to convert one set of a, 0 values represent- 
ing a critical point into another is needed. This is denoted as 
the Set Work. At the lower level, the work required to con- 
vert one a, 0 value in the set into another is needed. This will 
be defined as the Elemental Work. The Set Work is therefore 
the minimum amount of Elemental Work required to convert 
one set into another. Therefore, EMD can be used on the Set 
Work where the distance function is the Elemental Work. 


Definition 3 (Set Work) For two vector fields with feature 
distributions 


6 Feature Comparison via EMD 

A flow field can now be described by a set of a, 0 val- 
ues. To compare two flows, one approach is to find the clos- 
est match between the two sets of a, 0 values. The EMD 
algorithm provides this functionality. EMD is emphasized, 
since other techniques exist such as graph matching which 
takes into account connections between critical points. In 
the original description of EMD, terminology such as feature 
distribution, energy and work are used. The terminology is 
maintained for consistency and further information can be 
found in [14]. 

Earth Mover’s Distance computes the minimal amount of 
work required to transform one distribution to another. In 
the case of vector fields, the distribution can be represented 
as the set of ct, 0 values. 

Definition 1 (Feature Distribution) A feature distribution 
for a 3-D vector field is the set of sets of a and 0 values 


{{(nS 1 >/3} 1) ),(a( 1 ^< 1) ),(a«4 I) )}, 




and 


{{(^\0 1 {1) ),(a^ ) ,0 2 { \(a^\0 3 W )}, 


{(^ ) ,A w ),W,A w ),W,A w )}. 

. . . , {(a<"\ fB) ), (a<"\ $ n) )}} 

The amount of work necessary for transforming one vec- 
tor field into the other is defined as: 


( 2 ) 


( 2 ) (A( 2 ) 


Work set = X>MD e ( { (a^> , ) (a« , 0 ? ) , (4° , 0 ? ) } , 


i= 1 


7 Example: Lorenz Model 


EMD r=0 comparison for Lorenz Model o=1 0. 0=8/3, (KrS30 



Figure 3: EMD capturing topological changes to Lorenz 
Model 


where EMD e is the Earth Mover's Distance whose dis- 
tance function is the Elemental Work. 

Definition 4 (Elemental Work) For two vector fields with 
feature distributions 

and 

&),-■■ d<,P n )}- 

The amount of work necessary for transforming one 
vector field into the other is defined as: W ork e — 

v'Efci ((ai-a^ + Ofc-fl) 2 ). 

Notice that the definition for elemental work is identical 
to the work defined for a 2-D critical point defined in [1], 
Since we decompose a 3-D critical point into a set of 2-D 
critical points, EMD at this lower level is the same as the 
2-D case. To summarize, we can find the distance between 
two 3-D vector fields by extracting all 3-D critical points 
in the fields and representing them by a set of a, /3 values. 
Earth Mover’s Distance is used to find the minimum energy 
between the critical points just as in the 2-D case. However, 
the distance metric is defined as the EMD over the set of 
a, /? values whose distance metric is further defined as the 
elemental work. The elemental work is merely the Euclidean 
distance in a, 0 space and therefore is a metric. It can be 
shown that the Set Work is also a metric since EMD is a 
metric [14]. 


One application for field comparisons is in the area of me- 
teorology. Weather patterns can be searched for in a database 
to understand the development of a flow. E.N. Lorenz at- 
tempted to model thermally induced fluid convection in the 
atmosphere using the Navier-Stokes equations [15]. Us- 
ing two-dimensional motion, fluid heated from below and 
cooled from above under the effects of gravity produce cir- 
culation or convection rolls. This phenomenon is summa- 
rized in Equation 8. 

x 1 = o(y - x) (8) 

y‘ = rx — y — xz 

z' — —bz + xy 

x(t) represents a measure of the fluid velocity (amplitude of 
the convection motion), and y(t), z(t) represent measures of 
the spatial temperature distribution. The equations are in 
non-dimensional form where a is the Prandtl number (ratio 
of kinematic viscosity to thermal conductivity), b is a geo- 
metric factor and r is the Rayleigh number and is propor- 
tional to the temperature difference between the upper and 
lower regions of the system. These equations were one of 
the first to demonstrate chaotic behavior, and have resulted 
in over a hundred papers [16]. The resulting phase portrait of 
the system is in three-dimensions and is composed of simple 
critical points. 

In the original study by Lorenz er and b were fixed, and 
the Rayleigh number (i.e. the temperature difference be- 
tween the plates) was varied. When r is below 1 only one 
critical point exists. As the temperature difference increases, 
three critical points form and eventually for large enough r 
the entire system becomes unstable. Using a, /? space, the 
behavior of the system can clearly be seen and in fact the 
transformation from a stable system to an unstable system is 
continuous in this space. 

Beginning with r = 0, Equation 8 has one critical point 
at the origin. As seen in Figure 4a, the critical point is type 
Class I whose phase portraits are three attracting (stable) 
nodes. As r is increased to 1, the angle measured from the 
positive a axis reduces and therefore the nodes become less 
stable as it nears the saddle point The arrow in Figure 4a 
represents the critical point’s progression with increasing r. 
At r = 1, the angle is 135° and a degenerate node forms 
(Figure 4b). Increasing r further, the origin’s phase por- 
trait becomes a set of two saddles, and an attracting node 
(Figure 4c). As r increases, the saddle points approach the 
repelling node. 

In addition, for r > 1, two oth er critical points come 
into existence at ( ±y/b(r — 1), ±y/b(r — 1), r — 1). For 
r near 1, another set of stable attracting nodes form (Fig- 
ure 4d). However as r approaches 1 .346, the angle increases 


approaching an attracting star. Since the system is non- 
linear, the slightest perturbation causes the attracting star to 
become an attracting focus. Hence for 1.346 < r < 24.74 
(Figure 4e) an attracting focus forms. As r is increased 
to 24.76, the phase portrait changes gradually to a center. 
Hence, the observation of the unstable limit cycles forming 
around r = 13.926. Once r increases past 24.74, the an- 
gle increases to over 270° and a repelling focus comes into 
existence and the system is unstable. 

The change with temperature can also be understood by 
observing how the phase portrait for a particular r compares 
with the remaining phase portraits. Figure 3 plots the EMD 
values required to change the critical points for r = 0 into 
other critical points at other r values. As can been seen, 
for r < 1, a gradual increase in EMD occurs as the at- 
tracting node at the origin becomes a saddle. As soon as 
r > 1, two new critical points form causing a large jump in 
EMD. The jump is drastic since six additional (2-D) critical 
points (attracting nodes) must form. Once r > 1.346, the 
phase portrait changes from three attracting nodes, to one 
attracting focus. Hence, less energy is required to create two 
additional focii than six critical points. From r > 1.346, 
the EMD value slowly increases as the system continuously 
moves further from a stable system to an unstable one. Not 
only is the attracting focii becoming a repelling focii, but the 
saddles at the origin are approaching degenerate repelling 
nodes (Figure 4c). 


8 Conclusion and Future Work 

We have extended the feature based comparison method 
to three-dimensional vector fields. We have shown that the 
extension can be straight-forward if we use the property that 
a 3-D critical point can be decomposed into a set of 2-D 
critical points with planar phase portraits. In addition, the 
redefined distance function for EMD remains a metric. 

As with the 2-D case, connections between the critical 
points are not considered. Our experience has shown that for 
many cases this is sufficient, however, to reduce the num- 
ber of false positives and to provide a better distance be- 
tween two fields the connections should be taken into ac- 
count. We are investigating this aspect along with generaliz- 
ing this method to tensor fields. 

This method has also demonstrated its usefulness in un- 
derstanding complicated phenomenon such as the Lorenz 
model. The evolution of the thermal convection can be cap- 
tured with EMD. Since the system is represented by quanti- 
ties, fast searches can be easily constructed to locate partic- 
ular patterns. 
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Class Canonical Form Notes 


Phase Portraits 


IN PI 


1 


Xi 0 0 

0 A 2 0 

0 0 a 3 


Eigenvalues are real and 
distinct. Eigenvectors 
are linearly independent. 


A-Node, A-Node, A-Node 
R-Node, R-Node, R-Node 
A-Node, Saddle, Saddle 
R-Node, Saddle, Saddle 


Ql 

Q 2 
Q 3 


= 

\N^I) 

Ai +Aa 


V^2(Af+A§) 

A 2 +A 3 




01 = 

02 = 

03 = 


A 1 ~ A 2 

%/ 2 (A? +A 2 ) 

Ai -A3 

^/2(Af +A|) 
A2-A3 
x/2(A|+A|) 


Ai 0 O' 

0 A 2 0 

0 0 0 


2-D case 


A-Node 

R-Node 

Saddle 


Qi 


Ai+A? 
\/ 2 (Aj+Al) 
Q2 = 0 
Q 3 = 0 


01 = 


Ajr— Ag 


& = o 

& =0 


2 


Ai 0 0 

0 Ai 0 
0 0 A 2 


Eigenvectors are inde- 
pendent. 


A-Star, A-Node, A-Node 
R-Star, R-Node, R-Node 
R-Star, Saddle, Saddle 
A-Star, Saddle, Saddle 


Q2 

Q3 


Ql — ±1 

= A l+^2 

V^f+Af) 
— A 1 +A 2 

\/2(Ai-|-A|) 


0i=O 
Bo = 

\Z 2 (Ai+a!) 
/L, — ^1-^2 

^ \/2(A?+Al) 


Ai 0 O' 

0 Ai 0 

0 0 0 


2-D case 


A-Star 

R-Star 


Qi = ±1 0i = 0 

q 2 = 0 02 = 0 

q 3 = 0 03 = 0 


0 0 0 
0 0 0 
0 0 A 2 


1-D case 


ignore 


3 


Ai 1 0 

0 Ai 0 
0 0 A 2 


Eigenvalues are real. 
1 linearly dependent 
eigenvector. 


A-Log Star, A-Node 
R-Log Star, R-Node 
R-Log Star, Saddle 
A-Log Star, Saddle 


Ql = ±1 


Q3 


q 2 = 0 

~ v/2(Af +A1) 


01=0 
02=0 
03 = 

%/2(Af+Al) 


Ai 1 O' 

0 Ai 0 

0 0 0 


2-D case 


A-Log Star 
R-Log Star 


Ql — ±1 0i — 0 

0.2 = 0 02 =0 

03 = 0 03 =0 


0 1 0 
0 0 0 
0 0 A 2 


1-D case 


ignore 


Table 2: Phase Portraits for Classes 1,2,3. Legend: A- : Attracting R- : Repelling Ai \ 2 ^ A 3 ^ 0, Ai > A 2 > A 3 


Class Canonical Form Notes 


Phase Portraits 


ll Q 


11/311 


4 


Ai 

0 

0 


0 

Ai 

0 


0 

0 

Ai 


Eigenvalues are real. 
Eigenvectors are linearly 
independent. 


A-Star, A-Star, A-Star 
R-Star R-Star, R-Star 


Qi = ±1 / 3 i = 0 

0-2 = ±1 /?2 = 0 

Q3 = ±1 /S3 = 0 


5 


Ai 1 0 

0 Ai 0 
0 0 Ai 


Eigenvalues are real. 
One linearly dependent 
eigenvector. 


A-Log Star, A-Star 
R-Log Star, R-Star 


Ct! = ±1 Pi = 0 

Ct 2 = 0 P2 = 0 

«3 = ±1 Pz = 0 


6 


Ar 1 0 

0 Ar 1 
0 0 Ai 


Eigenvalues are real. 
Two linearly dependent 
eigenvectors. 


A-Log Star 
R-Log Star 


Ql — ±1 Pi = 0 

Q2 = 0 Pi = 0 

Q 3 = 0 P 3 = 0 


7 


a —P 0 
Pa 0 
0 0 A 3 


A-Focus 

Two complex and one R-Focus 
real eigenvalues. Center 


- 0 a. — 0 

\J a. 2 ~- 8 2 a 2 -rff 1 

a 2 = 0 P 2 = 0 

«3 = 0 p 3 = 0 


Table 3: Phase Portraits for Classes 4-7. Legend: A- : Attracting R- : Repelling Ai ^ A 2 ^ A 3 ^ 0, Ai > A 2 > A 3 
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Figure 4: Lorenz Model for Fluid Convection depicted in a, 3 space, a — 10, b = 8/3. AN:Attracting Node, 
DAX:Degenerate Attracting Node, SiSaddle, AF:Attracting Focus. RF:Repelling Focus 
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1 Abstract 

The previous work [1] introduced a quantitative method 
for comparing 2-D vector fields based on the number and 
type of critical points that comprised the field. However, 
the arrangement of the critical points was ignored poten- 
tially causing two very different fields with the same type of 
critical points but different streamline connections to be de- 
tected as similar fields. This paper improves the comparison 
technique by considering the connections between critical 
points, thereby improving the representation of the vector 
field. 


2 Introduction 

Comparing vector fields is an important analysis tech- 
nique used in the sciences. The comparison of fields aids 
in understanding and can improve mathematical models of 
phenomena. For example, a database containing historical 
weather patterns can be matched with the present day’s pat- 
tern to better understand the evolution of the flow. In fluid 
mechanics, simulated flow fields can be validated with ex- 
perimentally acquired data sets or with other models. Unfor- 
tunately, many of the techniques rely on rendering to a com- 
puter image and visually comparing the fields. This does not 
lend itself easily to automated comparison techniques, and 
certain fields may deceive the human eye. Data comparison 
techniques operate directly on the raw data but requires a 
common domain, grid alignment, and interpolation [2]. Fea- 
ture based comparisons are often problem specific, i.e. in 
fluid mechanics features include shock surfaces and vortex 
cores. These features are typically qualitatively compared. 

Recently in [1], a method for quantitatively comparing 
vector fields based on topological features was introduced. 
A vector field was characterized by the number and type of 
critical points. The pattern formed about the critical point 
was uniquely represented by an ( a,/3 ) coordinate. One of 


the properties of the (a, 0) space was that Euclidean dis- 
tances were valid and could be used as a difference measure 
between two critical points. Intuitively, the ( a,/? ) space 
describes the change of 2-D flow patterns as a progression 
around the unit circle. Traditionally, phase-portraits are clas- 
sified on p — q charts [3], Though other mappings are pos- 
sible, the (a, 0) representation provides a logical grouping 
and progression from one critical point to another thereby 
being attractive for comparing critical points. For exam- 
ple, vortical behavior is grouped together (/ 3 < 0) in the 
lower half of the unit circle. From dynamics it is known 
that a center pattern (0 = — 1) is structurally unstable (non- 
hyperbolic) and therefore is represented as a single point 
in this space. Any slight perturbation (typically introduced 
from non-linearities in the system) can change the pattern to 
an attracting or repelling focus (0 ss 1). Therefore the dis- 
tance between these focii can be small as found for nearly 
conservative systems. The Earth mover’s distance (EMD) 
was used as the matching algorithm to minimize the cost of 
work between the fields. The work in this case was merely 
the Euclidean distance between two (a, 0) points. However, 
the connectivity between critical points was not used. This 
could be problematic since two vector fields with similar 
critical points but significantly different connections led to 
similar EMD values and therefore were false positives. For 
a topologically accurate representation, it is necessary to in- 
clude not only the number and type of critical points but the 
set of connecting streamlines. 

Two methods are presented that define metrics for com- 
paring the number, type, and connections between vector 
fields. The first technique is based on graph theory where 
the vector field is represented using an attributed, relational 
graph (ARG) [4], The second method is based on EMD with 
a redefined feature distribution that captures the graphical 
structure of the vector field. The trade-off between the meth- 
ods is accuracy versus speed. The paper is concluded with 
an application demonstrating the advantages over the origi- 
nal method. 


3 Representation of Vector Fields using 
Attributed, Relational Graphs 


nents: 


An attributed relational graph is composed of a set of 
property based nodes and the relations between the nodes 
represented by edges or connections. In the case of vector 
field topology, the nodes are the critical points with ( a,0 ) 
as the characteristic property and the connecting streamlines 
are the edges that relate the critical points. To compare 
two attributed graphs, G, and Gj, a similarity or distance 
measure, D(Gj,Gj) is defined between them. The measure 
which can be thought of as a cost to transform G ,• to Gj 
involves both the cost of transforming the properties of the 
nodes and the relations of the graph. A desirable property 
for D(Gi, Gj) is that it be a metric [4]. 

The process of transforming the properties of the node 
(node-matching) can be interpreted as the total cost, Cn, re- 
quired to transform a set of nodes {p\ of G t to the set 
of nodes {<?} of Gj. Therefore, given a particular node- 
pairing, x , where f n (pi,Qj ) is the cost of transforming a 
pair of nodes (pi,qj), the total cost can be written as Cn = 

fniPiilj) where the summation is over all node pairs. 
For vector fields, the cost of node matching is merely the 
Euclidean distance between (a, 0) points [1]. If two vector 
fields contain an unequal number of critical points, normal 
points corresponding to ( a,0 ) = (0,0) are added. Normal 
points correspond to node insertion and deletion and have a 
uniform cost of unity since the distance from a normal point 
to any point on the unit circle is 1 . 

The relation between two graphs is equivalent to graph 
isomorphism which quantifies how close one graph struc- 
ture is to another. In order to equate two graphs, node inser- 
tion and deletion as well as edge insertion or deletion are re- 
quired. The edge relations between critical points for a vec- 
tor field are obtained by numerically integrating the stream- 
lines connecting the critical points [5). The edges can then 
be related by an adjacency matrix, A,j [4], This nxn matrix 
where n is the number of nodes (critical points) is assigned 
a value of 1 if critical point i is connected to critical point j, 
otherwise a 0 is assigned. Note that the adjacency matrix can 
represent only one connection between two nodes and has 
implications to be discussed. To compare graphs Gj and Gj, 
the corresponding adjacency matrix is compared component 
wise. The node-matching procedure described above guar- 
antees that the size of all adjacency matrices are the same 
since normal points are added to the flow field. Since the 
graph is undirected (i.e. node i is connected to node j, and 
vice versa), the adjacency matrix is symmetric and only half 
the matrix need be considered). The overall cost involved to 
convert vector field Gj with a configuration x and adjacency 
matrix A to another vector field, Gj with matrix A! is the 
sum of the relational-matching and node-matching compo- 


n n 

Ds(x) = Y,Y.\ A n- A 'ij\ 

i=l j=i 

+ ^y(Q i -a J ) 2 + (A-/3 J ) 2 (1) 

And the distance measure is defined as the minimum over all 
possible node pairings and the corresponding permutation of 
the adjacency matrix, 

D(Gi,Gj)= m i n D s (x) (2) 

3.1 Comments on the ARG Technique 

Since all permutations of the adjacency matrix need to be 
calculated, the overall complexity is 0(n! n 2 ) [5] which is 
unacceptable for large problems. However, significant re- 
search in the area of graph matching exists and is currently 
being investigated [6]. The adjacency matrix described is 
binary and does not take into account graphs with parallel 
edges. In the case of a saddle point, 2 streamlines may con- 
nect to a single critical point However, non-binary matrices 
can be used and is being investigated as well [7]. 


4 Extending EMD to Graph Matching 

As forementioned, an immediate disadvantage to the cur- 
rent implementation of ARGs is the time required for com- 
parison and therefore is an area of current investigation. 
However, a readily extendible solution is to modify the sig- 
nature defined in [1] to account for the connections between 
critical points. The EMD method in the worst case requires 
2” iterations, but in practice has been shown to be solved in 
polynomial time (cubic) [8]. 

In the previous work, a feature distribution for a vector 
field is the set of (a, 0) values associated with the critical 
points [1]. A new feature distribution is defined to be a set 
of (a, 0) values and the (a, 0) values of the immediate set of 
connections. This approximates a graph structure with n tree 
structures. The decomposition is not unique and therefore 
cases can be contrived such that two different graphs can 
be decomposed into the same tree structures [5]. However, 
given the relative position and types of critical points, there 
are a limited number of ways that streamlines can join the 
points, and as will be shown in the application section, in 
practice this method is comparable to ARG. 

The work required to convert one vector field into another 
is defined to be the Set Work. 

Definition 1 (Set Work) For two vector fields with feature 


distribution 


{{(a 1 ,/3i),(a ( 1 1) ,/3 1 (1) ),...,(Q ( 1 ni) ,/3i ni) )}, 

. . . {(a n ,0 n ), (c&Kffl), ■ ■ • , (a£*"\ /&"">)} } 

and 

{(Q 2 ,d 2 ),(4 1) ,^ 1) ),...,(4 n2) ,^ n2) )}, 

. . . {(d„ J n , (fiW W), . . . , (a<""> , /?("">)} } 



the amount of work necessary for transforming one vector 
field into another is defined as: 

n 

Workset = Work E s(,{(pcu0i)}, {(d»,&)}) 

i= 1 

+ min(Work B s){{((Xi 1) ,0i 1) ),--- , (a^ Bi) , / S|" i) )}, ... , 

where WcrrkES is the Elemental Set Work. 

The min( W or kgs) is the minimum amount of elemental 
work. This can be found using EMD. The Elemental Set 
Work remains the same and is defined in [1]. Therefore, 
the distance measure is defined as the minimum of the set 
work and is computed using EMD. To summarize, EMD is 
run at two levels. The first is at a global level where EMD 
minimizes the set work required to convert one vector field 
into another. The set work which converts one critical point 
into another requires the use of EMD to minimize the work 
required to convert one critical point’s children into another 
critical points’s children. 

4.1 Comments on Extending 
EMD Method 


Figure 1 : Attributed, Relational Graph compared with Orig- 
inal EMD method. Note differences in frames 3-7,13 and 
41. 
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Figure 2: Attributed, Relational Graph compared with Ex- 
tended EMD method. Note that the trends agree. 


As mentioned previously, the representation of a graph 
through a tree structure is not necessarily unique and there- 
fore can expect false positives. However, it may be possi- 
ble to determine up front whether a particular graph has a 
unique tree structure representation. In this case, the method 
works nearly as well as ARG. The results will not be the 
same, since a tree structure does not distinguish clearly be- 
tween the cost of node-matching and graph-matching as it 
does with ARG. In the limiting case, when no connections 
are present, the extended method produces the same results 
as the original EMD algorithm. Therefore, a lower limit can 
be established. The extended method is a metric since the 
Set Work is a composition of two metrics. Rubner has proven 
that EMD is a metric if the distance measure is a metric. [8]. 




Figure 3: Graph comparison of two vector fields. Region 1 
in Frame 4 extends to the upstream region of the airfoil. 


5 Applications 

To demonstrate the effectiveness of ARG and extended 
EMD, a comparison is made to the original EMD algorithm 
on the test case of a flow past a 2-D airfoil at —90° an- 
gle of attack. Rogers and Kwak computed this flow us- 
ing an incompressible, time accurate Navier-Stokes code at 
Re=200 [9]. Fifty time steps were computed capturing sev- 
eral periodic vortex shedding cycles. Figure 1 plots the EMD 
of frame 1 versus the remaining 49 frames. The plots have 
been normalized by the largest EMD values so the tech- 
niques may be compared. Intuitively, Figure 1 represents 
how much vector field 1 differs from the remaining 49 fields. 
Both methods capture the periodicity within the flow indi- 
cated by the valleys-every 16 to 17 frames. However, the 
ARG algorithm differs from the original EMD in frames 
3-7. From the topological illustrations in Figure 3, we see 
that contrary to EMD’s computation frames 1 and 4 are not 
equivalent. Region 1 in frame 4 extends to the upstream 
region of the airfoil. Also a new region is present start- 
ing in frame 13. Note that the original EMD method used 
only simple critical points and did not include the wall at- 
tachment/detachment points. In order to preserve the graph- 
ical structure, the separation points have been represented 
by normal points (a, 0 ) = (0, 0), but will be properly clas- 
sified in future papers. The number and type of attach- 
ment/detachment points do not change with time (always 2 
attachment and 2 detachment points) and so the methods are 
comparable. The cycle is consistently completed every 16 to 
17 frames. However, one marked difference occurs in frame 
41 which has an energy level lower than its counter parts in 
frames 8 and 25. This is due to frame 41 missing a sad- 
dle/repelling focus pair, and therefore the graph structure is 
closer to frame 1 . Hence the graph method has indicated a 
possible problematic area in the solver. 

Finally, Figure 2 compares the ARG method with the ex- 
tended EMD method. Immediately, it can be seen that both 
methods are in agreement with one another. The amplitudes 
are not expected to agree since each algorithm weighs the 
significance of the nodes differently. 
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6 Conclusions 

We have demonstrated the effectiveness of including con- 
nectivity information in the representation of vector field 
comparisons. By further considering the connecting fines 
forming the separatrices of the vector field, we are able to 
improve the quantitative measure between fields. For the 
airfoil data, we have shown that both the attributed graph 
method, and the extended EMD method can be used as an 
effective investigative and diagnostic tool. The extended 
EMD method is considerably faster than the brute force 
ARG method and in practice should be considered as an ef- 
fective first pass tool. 
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1 Abstract 

We study the topology of 3-D symmetric tensor fields. 
The goal is to represent their complex structure by a simple 
set of carefully chosen points and lines analogous to vector 
field topology. The basic constituents of tensor topology 
are the degenerate points, or points where eigenvalues are 
equal to each other. First, we introduce a new method for 
locating 3-D degenerate points. We then extract the topo- 
logical skeletons of the eigenvector fields and use them for a 
compact, comprehensive description of the tensor field. Fi- 
nally, we demonstrate the use of tensor field topology for 
the interpretation of the two-force Boussinesq problem. 

2 Introduction 

Second-order tensor fields have applications in many 
areas of physics, such as general relativity, fluid flows 
and mechanical properties of solids. The wealth of mul- 
tivariate information in tensor fields makes them more 
complex and abstract than scalar and vector fields. Vi- 
sualization provides a means to gain new insights from 
these rich data sets. 

The most natural way to visualize a symmetric 3-D 
tensor field is through its eigensystem, i.e., eigenvalues 
and eigenvectors. A continuous representation of the 
tensor field is obtained by tracing the trajectories of its 
eigenvectors. These trajectories are called hyperstream- 
lines [1,2]. The difficulty with such an approach is how 
to capture the structure of the 3-D domain while lim- 
iting the number of hyperstreamlines to a minimum in 
order to avoid visual clutter. The problem can be sig- 
nificantly simplified by taking a topological approach, 
similar to the one used in visualizing vector fields [3]. 
Degenerate points, defined as points where eigenval- 
ues are equal to each other, are the basic singularities 
underlying the topology of tensor fields. Eigenvectors 
never cross each other except at degenerate points. In 


the past, research has been conducted in the area of 
two-dimensional tensor fields [1, 2]. We live, however, 
in a three-dimensional world, and therefore, it is impor- 
tant for us to understand the underlying physics of this 
world. In this paper, we describe a new method for lo- 
cating degenerate points along with the conditions for 
classifying them in three-dimensional space. We also 
discuss some topological features of three-dimensional 
tensor fields, and interpret topological patterns in terms 
of physical properties. 

3 Theoretical Background 
3.1 Definitions 

Definition 1 (Second-Order Tensor Field) Let 

C(X, Y) be the set of all the linear transformations of 
the vector space X into the vector space Y, and let E be 
an open subset o/R n . A second-order tensor field T(x) 
defined across E is a mapping T : E —> £(R m ,R ra ) 
that associates to each vector x of E a linear trans- 
formation of R m into itself. If R m is referenced by a 
Cartesian coordinate system, T(x) can be represented 
by m? Cartesian components Tij(x), i,j = l,...,m, 
that transform according to 

m 

Tij = ^2 PipPjqTpq ( 1 ) 

p,q= i 

under an orthonormal transformation /? = {0ij} of the 
coordinate axes. [2] 

In a Cartesian coordinate system, a 3-D tensor field 
takes the following form: 

/ T u (x,y,z) T X2 (x,y,z) T 13 {x,y,z) \ 

T (x) = I T 2x (x,y,z) T 22 (x, y, z) T 23 (x,y,z) I 
V T 31 (x,y,z) T 32 (x,y, z) T 33 (x,y,z) ) 


( 2 ) 


Definition 2 (Hyper streamline) A geometric primitive 
of finite size sweeps along the longitudinal eigenvector 
field, vi, while stretching in the transverse plane under 
the combined action of the two transverse eigenvectors, 
v t i and v t 2 - Hyperstreamlines are surfaces that envelop 
the stretched primitives along the trajectories. We refer 
to hyperstreamlines as “ major ”, “ medium ” or “minor” 
depending on the corresponding longitudinal eigenvec- 
tor field that defines their trajectories and color hyper- 
streamlines by means of a user-defined function of the 
three eigenvalues, usually the amplitude of the longitu- 
dinal eigenvalue, [f] 

Definition 3 A degenerate point of a tensor field T : 
£'—>■£ (R m , R m ), where E is an open subset o/R m , is 
a point xq E E where at least two of the m eigenvalues 
of T are equal to each other. [4] 

3.2 Locating Degenerate Points 


A three-dimensional symmetric tensor field (Equa- 
tion (2)) has 6 independent variables, therefore various 
types of degenerate points may exist. These types cor- 
respond to the following conditions: 


Ai (T 0 ) = A 2 (T 0 ) > A 3 (T 0 ) (3) 

Ai (To) > A 2 (To) = A 3 (T 0 ) (4) 

Ai (T 0 ) = A 2 (T 0 ) = A 3 (T 0 ) (5) 


The characteristic equation of a 3-D symmetric tensor 
can be expressed in the following form 


A (A) 


Tu — A Xi 2 Ti 3 
Ti 2 T 22 — A r 23 
Ti 3 T 23 T 33 — A 

—A 3 T czA 2 T bX T c 


( 6 ) 


where a, b and c are composed of the 6 independent 
tensor components. The condition for the existence of 
a degenerate point is that both A (A (To)) and its deriva- 
tive are zero. 


A (A (;co)) — — A 3 + aX~ -t- 6A + c — 0 . . 

^(A(z 0 )) _ _ 3X 2 + 2a\ + b- f ' 


0 


As a result, we obtain the following conditions cor- 
responding to Equations (3, 4, and 5) respectively: 


0 . . 2a 3 + 9a6 + 2d 3 / 2 

Bi(x,y,z) = — bc = 0 (8) 


From the expressions for B\, S 2 and B 3 , we de- 
termine that: B\ ( x,y,z ) = 0 is a maximum for B i, 
f? 2 ( x , y, z) — 0 is a minimumfor f? 2 and B 3 ( x , y, z) = 0 
is a maximum for B 3 . 

Now the problem is to find extrema in a 3D contin- 
uous field from the discrete experimental data sets. On 
a 3-D discrete mesh, the search for the various extrema 
is conducted by processing one grid cell at a time for 
each spatial function. 

This method can successfully locate the points of 
triple degeneracy. It is especially useful when extended 
to locate points of double degeneracy where the local 
tensor appears in the diagonal form only when trans- 
formed into its eigenvector space. 

3.3 Separating Surfaces 

For second order tensor fields, in most cases, the 
eigenvector fields in the vicinity of a degenerate point 
can be described in terms of three types of angular sec- 
tors: hyperbolic, parabolic and elliptic sectors. It can 
be proved that in a 2-D tensor field, at a simple de- 
generate point, there are only one or three hyperbolic 
sectors, and no elliptic sectors [2], Correspondingly, we 
call the degenerate point a wedge point when it has only 
one hyperbolic sector and possibly one parabolic sector 
or a trisector point when it has three hyperbolic sectors 
[ 2 ]- 

The classification of degenerate points in 2-D tensor 
fields [2, 5] can then be extended to 3-D tensor fields. 
The building blocks are the fundamental elements as 
defined for 2-D [2], However, the separating surfaces in 
3-D tensor fields have a general structure as they could 
appear at various angles. Each of the surfaces is char- 
acterized by patterns similar to those of hyperbolic or 
parabolic sectors and is bounded by hyperstreamlines 
that are emanating from or terminated at the degener- 
ate point. Consequently, a point of triple degeneracy 
can be classified by the number and type of separating 
surfaces surrounding it. 

In Figure 1 we show the eigenvector patterns in the 
vicinity of a point of triple degeneracy with 4 bound- 
ing hyperstreamlines. These hyperstreamlines form 6 
hyperbolic separating surfaces. Figure 2 shows a point 
of triple degeneracy with only 3 bounding hyperstream- 
lines which form 2 hyperbolic separating surfaces and 
one parabolic surface. 

The trajectories on the surfaces are locally 2-D, while 
off the surfaces they are fully 3-D and are determined 
by their closest surface. 


b 2 (x, y, z ) = 


2a 3 + 9ab - 2 d 3 / 2 
27 


= 0 


B 3 (x, y, z) = a 2 + 36 = 0 


4 Topology of 3-D Tensor Fields 

We choose the elastic stress tensor induced by 
(10) two compressive forces on the top of a semi- infinite 



Figure 1: A point of triple degeneracy with 6 hyperbolic 
separating surfaces. 



Figure 2: A point of triple degeneracy with 2 hyperbolic 
separating surfaces and one parabolic surface. 

plane [6] to illustrate the advantages of using topo- 
logical skeletons in visualizing 3-D tensor fields. In 
principle, hyperstreamline trajectories of the stress ten- 
sor show the transmission of forces inside the mate- 
rial. Figure 3 shows two hyperstreamlines correspond- 
ing to the most compressive eigen-direction, the mi- 
nor eigenvector F3. The two forces, indicated by the 
arrows, act on the surface at Pi = (0.5, 0.0, — 1.05) 
and P2 = (—0.5, 0.0, —1.05) in the +z direction (down- 
ward). The domain of interest (described by the bound- 
ing frame) extends between (—1.0,-1.0,-1.114367) 
and (1.0, 1.0, 0.0) so it includes the key features of the 
stress tensor field, i.e. , the degenerate points. It is as- 
sumed that the region where z < —1.05 is in tension 
and that no stresses are transferred across the plane 
z — —1.05. The color of the hyperstreamlines encodes 
the magnitude of the most compressive eigenvalue, A3, 
while their cross section encodes the magnitude and 
direction of the transverse eigenvectors. The hyper- 
streamlines converge toward regions of high stresses 
where the forces are applied. Note the sharp change 
in color and cross section size of the hyperstreamlines 
as they approach the acting points of the forces. 

Analysis reveals that the tensor field contains two 
points of triple degeneracy and that these points re- 
side on the surface of the semi-infinite plane. More- 
over, the eigenvalues at these points ( the location 
of which is given by: = (0.0, 0.5, -1.05), D 2 = 

(0.0, —0.5, —1.05)) are equal to zero. This means that 



Figure 3: Stress tensor induced by two compressive 
forces; minor hyperstreamlines 



Figure 4: Stress tensor induced by two compressive 
forces; major hyperstreamlines 

these points are stress free, a fact that can be verified by 
an examination of the stress equations. We have there- 
fore acquired physical insight into the stress tensor field 
just by examining a basic topological feature, a point 
of triple degeneracy. 

Figure 4 shows hyperstreamlines that are obtained 
by tracing the major eigenvector field. The location 
and direction of the forces are indicated by the arrows 
and the location of the points of triple degeneracy are 
marked by spheres. The hyperstreamlines are presented 
with a constant cross section to avoid visual clutter 
resulting from the high eigenvalues in the vicinity of 
the points of the acting forces. They are, however, 
still color encoded by the major eigenvalue. Each of 
the 2 degenerate points has 4 bounding hyperstream- 
lines(separatrices), three of which lie on the surface 
z = —1.05 in a trisector pattern and the forth, which 
is pointing in the +z direction, connects the points of 
triple degeneracy, and delineates one of the two symme- 




5 Conclusions 



Figure 5: Stress tensor induced by two compressive 
forces; minor hyperstreamlines 



Figure 6 : Stress tensor induced by two compressive 
forces; medium hyperstreamlines 


try planes (the other goes through the points of action 
of the forces). 

To further clarify the tensor topology, the skeletons 
of the minor and medium hyperstreamlines are pre- 
sented in Figures 5 and 6 respectively. We can see 
from Figure 5 that the minor hyperstreamlines form a 
trisector-point like pattern in the vicinity of the points 
of triple degeneracy. They also indicate that a locus 
of points of double degeneracy (A 2 = A 3 ) connects the 
points of triple degeneracy. This is evident from the two 
trisector points that lie in the symmetry planes just be- 
low the points of triple degeneracy. The existence of the 
line of double degeneracy is further verified by noting 
the two points of double degeneracy in the skeleton of 
the medium hyperstreamlines (Figure 6 ). 


In this paper, we applied novel methods to determine 
the topology of tensor data sets, and made use of ad- 
vanced representations to determine the significance of 
degenerate points and topological skeletons in terms of 
physical features. 

By extracting the geometric structure of tensor data, 
we produce simple and austere depictions that allow 
observers to infer the behavior of any hyperstreamlines 
in the field. It enables important elements of 3-D stress 
distribution to be envisaged without visual clutter. 

Degenerate points represent the singularities of the 
tensor field. In the 3-D elastic stress tensor case we were 
able to identify points of zero stresses with triple de- 
generate points and to illustrate transmission of forces 
inside the material. 
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Abstract 

Combinatorial topology, also known as “rubber 
sheet geometry”, has extensive applications in ge- 
ometry and analysis, many of which result from con- 
nections with the theory of differential equations. A 
link between topology and differential equations is 
vector fields. Recent developments in scientific visu- 
alization have shown that vector fields also play an 
important role in the analysis of second-order tensor 
fields. A second-order tensor field can be transformed 

* Graduate Student, Physics Department 

t Graduate Student, Aeronautics and Astronautics 

Department 

* Professor, Electrical Engineering and Aeronautics and 

Astronautics Departments, Member AIAA. 

§ Lecturer, Faculty of Aerospace Engineering, Member 

AIAA. 

Copyright © 1997 by Authors. Published by the 
American Institute of Aeronautics and Astronautics, Inc. 
with permission. 


into its eigensystem, namely, eigenvalues and their 
associated eigenvectors without loss of information 
content. Eigenvectors behave in a similar fashion to 
ordinary vectors with even simpler topological struc- 
tures due to their sign indeterminacy. Incorporating 
information about eigenvectors and eigenvalues in a 
display technique known as hyperstreamlines 1 reveals 
the structure of a tensor field. To simplify an often 
complex tensor field and to capture its important fea- 
tures, the tensor is decomposed into an isotropic tensor 
and a deviator. A tensor field and its deviator share 
the same set of eigenvectors, and therefore they have 
a similar topological structure. A deviator determines 
the properties of a tensor field, while the isotropic 
part provides a uniform bias. Degenerate points are 
basic constituents of tensor fields. In 2-D tensor fields, 
there are only two types of degenerate points; while 
in 3-D, the degenerate points can be characterized in 


a Q’-R’ plane. Compressible and incompressible flows 

fields arising from gravitation and electromagnetism; 

share similar topological features due to the similarity 

the velocity vectors of a fluid motion, such as the 

of their deviators. In the case of the deformation 

atmosphere; and gradients, such as the pressure gra- 

tensor, the singularities of its deviator represent the 

dient on a weather map or the height gradient of a 

area of the vortex core in the field. In turbulent 

relief chart are all elements of vector fields. Recent 

flows, the similarities and differences of the topology 

developments in scientific visualization have shown 

of the deformation and the Reynolds stress tensors 

that vector fields and their topological structures also 

reveal that the basic eddie-viscosity assumptions have 

play a very important role in analyzing second-order 

their validity in turbulence modeling under certain 

tensor fields. 

conditions. 

Second-order tensors are fully represented by their 

Introduction 

eigenvectors and associated eigenvalues. 


ii 

H 

Combinatorial topology is a branch of geometry. It 

Where A, and e] (i — 1,2,3) are eigenvalues and 

studies the properties of figures that endure when the 

eigenvectors of the tensor T, respectively. The A,’s 

figures are subjected to continuous transformations. A 

represent all the amplitude information while the 

topological property of a figure is a property possessed 

e^’s represent all the directional information of T. 

alike by the figure and all its topological equivalents. 

Visualizing a tensor field is equivalent to visualiz- 

In combinatorial topology, complicated figures are 

ing its eigenvector fields. However, unlike a vector 

constructed from simple ones and their properties are 

field, Eigenvectors are vectors with sign indetermi- 

deduced from the simple figures. 2 

nacy. This remarkable feature distinguishes eigen- 

Combinatorial topology has vast applications to 

vector fields from ordinary vector fields and makes 

solving systems of differential equations. It may seem 

their topological features even simpler. The basic 

surprising that such superficially different subjects as 

constituents of tensor topology are degenerate points 

topology and differential equations could be related, 

where at least two of the eigenvalues are equal to each 

but research has shown that a link between these 

other. They play a role similar to critical points in 

two subjects is the concept of a vector field. Vector 

vector fields. As with vector fields there are numerous 

fields have many important applications. The force 

tensor quantities that are of interest to analyze and 
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visualize. In fluid mechanics, the Reynolds stress and 
the strain-rate (deformation) tensors are examples of 
tensors of high interest. This is in part due to the fact 
that eddie viscosity turbulence models are based on 
the assumption that the Reynolds stress tensor can 
be modeled as a function of the deformation tensor. 
A comparison of the topology of the two tensors 
would provide a means for validating the basic as- 
sumptions. Another aspect of the deformation tensor 
is the alignment of its medium eigenvector with the 
vorticity vector and scalar gradients. 3 4 It is believed 
that a stationary flow near the vortex core under goes 
rigid body motion implying that the components of 
the deformation tensor are zero. Computation and 
visualization of the phenomena is still an area of 
research and debate. Other tensors commonly used in 
the area of fluid mechanics are the velocity gradient, 
stress tensor, and the momentum-flux-density tensor. 

Theoretical Background 

Critical Points and Flow Field Topology 

A critical point is a point in the flow field where all 
three velocity components are zero and the streamline 
slope is indeterminate. The number and nature of crit- 
ical points in a vector field remain unchanged under 
a continuous transformation. These are topological 


properties of the system of differential equations. 

Historically, critical point theory has been used pri- 
marily to examine the solution trajectories of ordinary 
differential equations. (see Kaplan , 5 Pontryagin, 6 An- 
dronov,' and Minorsky 8 ) In a reference frame where 
a critical point is located at the origin, the motion in 
the flow is described by the leading terms of a Taylor 
expansion for the velocity field. The streamlines 
are, therefore, defined by the solution trajectories of 
three linear, coupled, first-order ordinary differential 
equations. The relationship between the properties of 
this 3x3 Jacobian matrix and the geometry of the 
solution trajectories is not trivial. 

A critical point can be classified according to the 
eigenvalues of its Jacobian matrix . 9,10 A positive or 
negative real part of an eigenvalue indicates a repelling 
or attracting nature, respectively. The imaginary part 
denotes circulation about the point. Accordingly, a 
critical point can be classified as an attracting node , a 
repelling node, an attracting focus, a repelling focus, a 
center or a saddle. 

For a complex flow field, the eigenvector planes are 
not always easily identified and located and it becomes 
necessary to use an analysis based on three invariants 
P,Q and R of a 3 x 3 Jacobian matrix. 11 All possible 
linear local flow trajectories of a moving continuum 
for both compressible and incompressible fluids can 
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be completely categorized in the space of P,Q and Ft. 
A set of surfaces can be defined in this space which 
defines boundaries between topologically distinct flow 
patterns and serves as a guide for identifying critical 
points. 11 

Tensor Field Analysis 

A tensor field in the real world is often very com- 
plex and difficult to analyze as a whole. Therefore, 
reducing the field to a simpler form is desirable. 

Decomposition of a Tensor Field 

Any given tensor can be decomposed into a sum of 
a deviator D and an isotropic tensor U : 

T = D + U (2) 

Definition 1 (Deviator) A tensor is a deviator D 
iff it is trace free, i.e. Trace (D) — 0, 

Definition 2 (Isotropic Tensor) 

A tensor is isotropic iff Uij = vhij . where v is a 
stretch factor. 

Definition 3 (Degenerate Point) A degenerate 
point of a tensor field T : E — * £(R m ,R m ), where E 
is an open subset of R m , is a point xq € E where at 
least two of the m eigenvalues of T are equal to each 
other. 


Definition 4 (Singular Point) A singular point in 
a tensor field is a point where all eigenvalues of a 
tensor vanish, in mathematical representation, it is a 
zero matrix. 

At any location in a tensor field, an isotropic 
tensor behaves the same in every direction, in other 
words, it is isotropic throughout the whole field. The 
topology of such a tensor is fairly simple and easy to 
deduce. 

A deviator, in contrast to an isotropic tensor, 
has a different behavior in all 3 principal directions 
except at a singular point where all of its components 
are zero. A general tensor and its deviator have 
the same set of eigenvectors; and degenerate points of 
a general tensor and singular points of its deviator 
occur at the same locations. 12 This means that the 
topology of a tensor field is identical to the topology 
of its de viator. A real tensor field is a deviator super- 
imposed onto an isotropic tensor. By subtracting 
the contribution of the isotropic tensor from the 
tensor field, the deviator becomes dominant. This 
allows a clear depiction of the topology and the fluc- 
tuations of the field without the disturbance of the 
sometimes dominant isotropic contribution. In the 
following classification of a 3-D tensor field topology, 
analysis of the deviator also provides a mathematical 
simplification that reduces the otherwise 3-D P-Q-R 
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space into a 2-D Q-R plane. 

2-D Tensor Fields 

Hyperstreamlines 

In the case of two-dimensional tensor fields, there 
are only two eigenvalues Ai and A 2 , and xo is a 

Similar to streamlines in a flow field, hyperstream- 

degenerate point iff Ai (xo) = A 2 (xo). 

lines are integral curves satisfying: 

Similar to critical points in vector fields, tensor 

f=*(*) <»> 

fields have different types of degenerate points that 
correspond to different patterns in their neighbor- 

Where s is a parameter measuring distance along the 

hoods. These patterns are determined by the tensor 

path. Unlike the velocity field, v represents one of the 

gradients at the degenerate points. 14-16 

eigenvectors of a tensor field T (x) obtained from: 

Consider the partial derivatives 

T(x 0 )F(x 0 ) = A(x 0 )t>(x 0 ) (4) 

- 1 d(Tn— T 22 ) 1 1 d(Tn — T 22 ) 

“ _ 2 dx u - 2 8y ... 

(5) 

at each location xo. 

O 

II 

3f 

II 

A hyperstreamline in 3-D is a trajectory that traces 

evaluated at the degenerate point Xo- In the vicinity 

along the longitudinal eigenvector field while stretch- 

of xo, the expansion of the tensor components to 

ing in the transverse plane under the combined action 

first-order is as follows: 

of the two transverse eigenvectors. We refer to hy- 

t il=T 22 ^ a ^ x frAy 

(6) 

perstreamlines as “major”, “medium” and “minor” 

r i2 « cAx + dAy 

depending on the corresponding longitudinal eigen- 

where (Ax, Ay) are small displacements from xo- 

vector field that defines its trajectory. 

An important quantity for the characterization of 

Topology of Tensor Fields 

degenerate points is the quantity 

The topology of a tensor field T(x) is the topology 

6 = ad — be (7) 

of its eigenvector fields t7,-(x). 13 Similar to critical 

The appeal of 6 arises from being invariant under 

points, degenerate points are the basic constituents 

rotation. When 6 < 0, the degenerate point has 

underlying the topology of tensor fields. Degeneracy 

three hyperbolic sectors. The pattern of eigenvector 

occurs at points where at least two eigenvalues are 

fields corresponding to the trisector point is shown in 

equal to each other. 

Figure 1. When 6 > 0,the degenerate point has one 
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In P-Q-R space, it can be shown that the solution 



trajectories only exist between surfaces Si and S 2 , 
which are, respectively, given by: 

2 P 3 + 9PQ + 2 (P 2 + 3 Q) 3/2 


27 


+ R = 0 (12) 


and 


2 P 3 + 9PQ -2 (P 2 + 3 Q) 


3/2 


27 


+ i? = 0 (13) 


Figure 1: Trisector (6 < 0) and wedge (6 > 0) points. 

hyperbolic sector. The local pattern corresponds to 
the wedge point represented in Figure 1. 

3-D Tensor Fields 

The eigenvalues of a 3-D tensor T can be obtained 
by investigating its characteristic equation. 


where: 


P = Tn + T22 + T3: 


Q = 


Tn 

Tl2 


Tn 

T13 


T 22 

T 2 3 



+ 



+ 



Tl2 

T22 


Tl3 

T33 


T23 

T33 


(9) 


( 10 ) 


Til T12 T13 

R = T\ 2 T22 r 23 (11) 

Tis T23 T33 

Similar to critical point analysis, one can see that 


On surface Si, Aj (xo) = A 2 (xo); and on surface 
S 2) A 2 (^ 0 ) = A 3 (*„). Triple degeneracy Ai (« 0 ) = 
A 2 (zq) = A 3 (x 0 ) occurs at points where Si and S 2 


• n p3 p 2 

meet, i.e. R — Q — 


A — Tn 

-Ti 2 

-T13 


1 

-T12 

A — T22 

— T23 

= A 3 + PA 2 + QA+R (8) 


-T13 

— T23 

a-t 33 




Since the tensor and its deviator have the same 
topology, it is sufficient to examine only the deviator 
D. By definition, D\\ + £>22 + P 33 = 0, and therefore, 
P — 0. The coefficient Q can also be presented as: 


(14) 


The characteristic equation now becomes: 
A 3 + Q'X + R' = 0 


(15) 


The P-Q-R space now reduces to the Q’-R’ plane. 
On this plane, the solution trajectories exist between 
curves L\ and £ 2 (shaded area). L\ and L 2 are given 
respectively by 

Hf) § 

and 


R!_ 

2 


Q'\ i 


the coefficients P, Q and R are all tensor invariants. 
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(17) 


Case Studies 



Figure 2: Q'-R' plane 


Double degeneracy occur on the curves L\ and Ln 
and triple degeneracy occurs where L\ and I 2 meet, 
i.e. <? = R' = 0. 

For double degeneracy cases, the pattern of the 
hyperstreamlines resembles that of a 2-D tensor 
field. For example, consider a degenerate point 
where At. ( xq ) = A 2 (To) > A3 (To). The tensor field 
is degenerate in the plane orthogonal to ^ 3 ( 10 ) within 
which locally two-dimensional patterns such as wedge 
points and trisectors can occur. 

However, the case of triple degeneracy is global 
and the structure of the eigenvector trajectories in 
the vicinity of the degenerate point is fully three 
dimensional. The general structure of the eigenvectors 
in the vicinity of a point of triple degeneracy is based 
upon the fundamental 2-D degeneracies in planes 
defined by the 3-D tensor and its expansion in the 
neighborhood of (To). 1, 


2-D flow past a hemisphere cylinder 

The concepts discussed above are illustrated by first 
visualizing the topology of the stress tensor of a 2-D 
flow past a cylinder. Fluid elements undergo compres- 
sive stresses while moving with the flow. Stresses are 
described mathematically by the stress tensor, which 
combine isotropic pressure and anisotropic viscous 
stresses (deviator). Both eigenvalues of the stress 
tensor are negative and the two orthogonal eigenvec- 
tors, vY and uT, are along the least and the most 
compressive directions, respectively. At a degenerate 
point, the viscous stresses vanish and both eigenvalues 
are equal to the pressure; degenerate points are points 
of pure pressure. In Figure 3, the texture' 1 represents 
the most compressive eigenvector of the stress tensor 
(rT). Color encodes the magnitude of the compressive 
force (A), from most compressive(red) to least com- 
pressive(blue). Red dots marked with W are wedge 
points; white dots marked with T are trisectors. " 

-The texture is created by a technique discussed in Refer- 
ences . 16 

“This part of the work has been done by Thierry Delmar- 
celle, for details refer to 16 



Figure 3: Stress in a flow past a cylinder; most Figure 4: Reversible momentum flux density tensor in 

compressive field. the flow past a hemisphere cylinder 


Vortical Flow about a Slender body of 
Revolution 


In the next example, hyperstreamlines are used 

to correlate several different physical quantities in 

fluid flow visualization. For the reversible part of 

the momentum-flux-density tensor 11^ = p6ij + pv{Vj 

( p is the mass density and v, and Vj are the velocity 

components in the i, j direction, respectively), one may 

correlate pressure p , velocity direction, and kinetic 

energy density |pt> 2 (p is the velocity magnitude). 

Indeed, the tensor field fF can be diagonalized as 

/ , \ 

p + pv 2 0 0 

n r = o P o 

0 Op 

The major eigenvalue of Il r is Ai — p + pv 2 and 
the associated unit eigenvector is aligned with the 
velocity vector. The other eigenvalues are degenerate 


V 


/ 


and equal to the pressure (A2 = A3 = p) in the 
whole space. It follows that only major eigenvectors 
can be traced. Their trajectories are tangent to the 
velocity field and correspond to streamlines of the 
velocity field. The tubes’ cross-sections are circular 
with a diameter proportional to the pressure p. 16 In 
Figure 4, color encodes the kinetic energy density. 
The direction of the incoming flow is 5° with respect 
to hemisphere axis, the Reynolds number is 14,000 
and the flow is incompressible. The detachment 
at the end of the cylinder is clearly visible. The 
pattern of hyperstreamlines indicates that momentum 
is transferred from the nose of the body to the end 
fairly uniformly with a globally decreasing kinetic 
energy, as shown by color variations. However, there is 
a sudden change of kinetic energy (color) and pressure 
(diameter) associated with a significant variation of 
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the direction of the first five tubes. 

As noted before, a tensor field and its deviator 
share the same topological structure. The importance 
of a deviator in tensor analysis can be demonstrated 

( , , n \ 
xy -x i + yf 0 

2r 

72 -x 2 + y 2 -xy 0 08) 

\ ° ° ° 

Here r is the distance from a point to the center of 

by the following discussion and illustration of the 

the vortex core and R is the radius of the vortex core. 

deformation tensor in a stationary flow field. 

It is virtually a 2-D tensor with major and minor 


eigenvalues having equal magnitude but opposite sign 

For an incompressible stationary flow the deforma- 

and the medium eigenvalue remains zero. The defor- 

tion tensor, defined as Def = has a zero 

mation tensor is discontinuous at r = R. The angles 

isotropic part and therefore is equal to its deviator. 

of separatrices can be calculated by using the tensor 

For a compressible flow, the deformation tensor is 

in the neighborhood of the vortex core Def (r > R). 16 

composed of a deviator superimposed on a non-zero 

However, there is no real solution for the angles. This 

isotropic tensor which represents the rate of expan- 

indicates that the major and minor eigenvector fields 

sion. Therefore, a deviator describes the topological 

are a pair of loci in the transverse plane while the 

structure for both incompressible and compressible 

medium eigenvector follows the direction of the vortex 

stationary flows. 

core. 


Figure 5 shows hyperstreamlines of the deformation 


Following the assumption that for a rotational flow, 

tensor of flow past a hemisphere-cylinder. The angle of 

inside the vortex core, the flow is purely rotational. 

attack is set at 19° with a Reynolds number of 445,000 

Assume that the flow advances in the z-direction and 

and a freestream Mach number of 1.2. There is a pair 

rotates about the z-axis while the velocity within 

of primary vortices in this flow. The eigenvector field 

the vortex core area is (—u>y,ux,v z ) (where u is the 

around each core forms a focus swirling opposite to 

angular velocity and v z is the axial velocity and is 

each other. Both minor and major eigenvector fields 

a constant for a stationary flow) and its deforma- 

form a ring pattern. The upper ring in Figure 5 is a 

tion tensor Def (r < R) becomes singular; outside 

minor hyperstreamline and the ring in the lower part 

the vortex core, the velocity is p- (— y, a:,0) and its 

which encloses the body is the major hyperstreamline. 

deformation tensor Def (r > R) is: 

The two hyperstreamlines along the body are the 
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Figure 5: Deformation tensor in a flow past a hemi- 
sphere cylinder at incidence 



medium eigenvectors and also define the direction of 
the vortex core. 


Figure 6: Deformation tensor in the near field of 
flow past a wingtip; major eigenvector field 


Near Field of a Wing-Tip Vortex 


Eddie viscosity turbulence models have been widely 
used in the simulation of turbulent flow, mainly for 
modeling attached flows. Under certain conditions, 
some of the models can be modified and extended 
to model separated flows around slender bodies of 
revolution at high angle of incidence. 

In the following example, the near field of a wing- 
tip vortex is examined to determine whether the basic 
assumptions behind the eddie viscosity turbulence 
models are valid. If not, the question as to why do 
they still assist in simulating the flow with relatively 
reasonable results needs to be answered. In an ex- 
periment that was conducted by Chow and Zilliac, a 



Figure 7: Deformation tensor in the near field of 
flow past a wingtip; minor eigenvector field 
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Figure 8: Reynolds stress tensor in the near field of a 
flow past a wingtip; major eigenvector field 



Figure 9: Reynolds stress tensor in the near field of a 
flow past a wingtip; minor eigenvector field 


wide array of measurements, including surface oil-flow 
visualization, laser-illuminated smoke visualization, 
surface pressure and velocity-field measurements by 
use of a 7-hole pressure-probe and triple-wire hot- 
wire anemometry, was completed for the flow over a 
rectangular wing with rounded tip, up to 0.67 chords 
downstream of the trailing edge. The angle of attack 
of the wing was set at 10° and a chord Reynolds 
number was chosen at 4.6 million based on the desire 
to study a fully turbulent vortex. ^ 

The primary objective of post measurement anal- 
ysis is to compare the deformation tensor with the 
Reynolds stress tensor. Topological methods are 
ideally suited for this purpose. Previously, methods 
for visualization of tensor fields ranged from com- 
paring the individual components of the tensors to 
comparing the eigenvalues and eigenvectors. However, 
these conventional methods have failed to determine 
whether the tensors are similar. Here we describe a 
new approach. 

The eigensystems are computed for both Reynolds 
stress and deformation tensors obtained from exper- 
imental data. Texture maps are used to display the 
fields in the transverse plane across the vortex core. 


' * for detailed information about this experiment, please 
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From the eddy-viscosity hypothesis, 

-pRij = P»T ^ (19) 

Where Rij is the Reynolds stress, ; s the mean 
deformation tensor and l>t is the kinematic eddy 
viscosity. The eigenvalues of the Reynolds stress 
tensor are equal to the negatively scaled eigenvalues of 
the mean deformation tensor. Therefore, the order of 
eigenvalues of two tensors are opposite to each other, 
i.e. Ain > A #2 > while Af/i < Af /2 < \uz- This 
leads to the fact that the major eigenvector field of the 
mean deformation tensor swirls in the same direction 
as the minor eigenvector field of the Reynolds stress 
tensor. Figure 6 and Figure 7 display major and 
minor eigenvector fields for the deformation tensor, 
respectively. Their patterns show a pair of loci which 
confirms the discussion above regarding the behavior 
of eigenvector fields of a deformation tensor in a 
vortex core area. Figure 8 and Figure 9 are major and 
minor eigenvector fields for the Reynolds stress tensor, 
respectively. Their patterns also show a pair of loci. 
However, the minor eigenvector swirls faster than that 
of the deformation tensor, and therefore has a tighter 
focus. Similarly, the major eigenvector swirls slower 
and has a looser focus than the minor eigenvector 
fields of the deformation tensor. The similarity of the 
eigenvector fields between the deformation tensor and 
the Reynolds stress shows that the basic assumption 


behind the eddie-viscosity turbulence model are valid 
to a limited extent because the fields of these two 
tensors have different swirling speeds. This indicates 
that the tensors are not aligned and therefore the 
model is too simple for capturing the complicated 
turbulent flow in the near field of the wingtip. It is 
suggested that the model can be improved by taking 
into account the variation in swirling speed. 
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Abstract 

This report discusses a method for extracting and visualizing shock waves from three dimen- 
sional computational data-sets. Issues concerning computation time, robustness to numerical 
perturbations, and noise introduction are considered and compared with other methods. Fi- 
nally, results using this method are discussed. 


Introduction 

The visualization of shock waves within transonic, supersonic, and hypersonic flows has been 
a persistent problem. The difficulty lies in accurately extracting the shock from a computed 
data-set and presenting a meaningful representation of the results. Accurately extracting and 
depicting these shocks can lead to improvements in aerodynamic design and an insight into 
the physics of the flow. Much work has been done on the extraction of features within large 
data-sets using various data visualization techniques such as isosurfaces, stream ribbons, 
and contour maps. Recently, work done by Pagendarm and Seitz [1] has been successfully 
applied towards the extraction of shock waves. Their technique relied on the fact that the 
density jump across a shock numerically smears across a narrow band. This can be seen in 
Figure la, which represents a 1-D example, where the £ axis is a location along the shock 
and p is the density. Since A£ is a rather small distance, the shock can be assumed to be 
located at the inflection point. This point can be located via the second derivative, |^f = 0. 
Figure lc depicts this derivative. Many numerical algorithms can be applied to extract this 
zero value. Pagendarm, and Seitz did not use the second derivative directly, since they were 
dealing with 3-D datasets; instead, they used u ■ V(u • Vp). There are a few difficulties 
with this method. From a numerical method standpoint, the noise is amplified for each 
derivative taken. Furthermore, this method requires the gradient to be taken twice, thereby 
producing many local minima and maxima which are detected by zero-search algorithms. 
These artifacts can be filtered with some computational effort. Another problem is that 
this method may cause false detection in free stream regions. The minute perturbations in 
numerical data in the free stream produce regions where the first and second derivatives 

'Professor, Electrical Engineering and Aeronautics and Astronautics Departments. 
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Figure 1: Density variation and its derivatives across a shock wave 

appear to cross zero, rather than remaining fixed at zero. This causes zero-search algorithms 
to tag these points. Pagendarm’s and Seitz’s algorithm, however, is general enough to be 
applied to a variety of problems concerning discontinuities. What we propose is a technique 
that takes advantage of shock attributes for a faster and cleaner extraction. Faster in that 
we minimize the operation count, and cleaner in that we reduce noise. 


Theoretical Background 


A shock represents a sudden change of fluid properties. Typically, a shock is witnessed when 
a body travels at supersonic speeds. The flow adjusts to a body by abruptly changing its 
pressure, density, and temperature. This abruptness is caused by the flow’s inability to 
sense the body. The Mach number M which is defined as the ratio of stream velocity to 
sonic velocity can be composed of two components namely one that is parallel M||, and the 
other that is perpendicular Mi to a shock wave. The fluid properties such as density, static 
pressure, and total temperature are distinct on either side of the shock and appear as a jump 
across the shock. We define Mi as follows: 


Mi 


Vi-Vi 


( 1 ) 


where a is the speed of sound and Vi is the velocity in the direction of the pressure gradient, 

VP 


V, 


V- 


|VP| 


( 2 ) 
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The transition of Mi across a shock progresses from supersonic, Mi > 1 to subsonic, Mi < 
1. Equation (1) assumes that the pressure gradient is normal to the shock. We use this 
assumption in our method for shock extraction. 

Several steps are required for the extraction of the shock. We first begin by computing 
Mi for the entire computational space using Equation (1). Mi is used, rather than Mach 
number, to ensure the capture of oblique shock waves. The computational space is composed 
of grid cells. Each cell, composed of eight neighboring grid points, is checked to ensure the 
region is in compression VP • V > 0. Shocks exists only in regions of compression; by 
removing expansion waves, further computation is not required. A further check to validate 
the existence of a shock within a cell, is to validate that the change in pressure ( 8P ) is within 
a safety factor of the Rankine-Hugoniot range. Regions where the pressure gradient changes 
direction abruptly, thereby producing a large angle between the velocity vectors, can cause 
Mi to drop below 1.0, and appear as a shock. Therefore, using the isentropic relation, 


9 * 'v 

6P = P* * M^ 

7 + 1 


7 ~ 1 
7+1 


( 3 ) 


where Mx is the largest in the computational set, we ensure that the pressure gradient is 
large enough for a shock to exist. If a cell meets these criteria, it is then passed through the 
marching cube algorithm [2]. Surfaces constructed from triangles of Mx = 1.0 are created. 
These triangles are further checked to ensure that the pressure gradient is in fact in the 
direction of the surface normal. This filters out regions near the body (in the boundary 
layer) that can be tagged as a shock since the velocity progresses from zero at the body 
(no slip condition) to free stream velocity, hence passing through Mx = 1.0. This method 
of shock detection requires only the first derivative VP to be computed, which introduces 
less noise and is computationally faster than other methods requiring a second derivative. 
Computation time is further reduced by making simple checks to ensure a region’s candidacy 
for a shock prior to extensive processing. This method is also useful for finding shocks in 
solutions that have not converged fully. Free stream perturbations will not be tagged since 
characteristics that are indigenous to flow fields are used. Therefore, this method is not ideal 
for general 3-D data-sets such as medical imaging data but works best for flow problems. 


Results 

Several test cases were run to ensure the validity of the method. Figure 2 depicts a lamda 
shock extracted from a data-set of an ONERA M6 wing 1 with an incoming Mach number, 
Moo = 0.84, angle of attack, a = 6°, and Reynolds number, Re = 761,000. As can be seen 
from the figure, we have a clean extraction of the shock with little to no artifacts. The 
extraction also reveals a region where the lamda shock is separated. This space informs the 
designer that the solution may not have converged fully, or a higher grid resolution is required 
in the area. As another verification of the shock algorithm, we extract shock features from 
a hemisphere cylinder at an angle of incidence of a = 19°, Moo = 1.2, and Re = 445,000). 
Figure 3 depicts the results. Two distinct shocks exist. The first is the bow shock which is 
upstream of the nose. The coarseness of the grid yields ’’holes” in the shock surface. Further 

1 ONERA (Office National D’etuedes Et De Recherches Aerospatiales) wing designed for studies of 3-D 
flows of low to transonic speeds at high Reynolds numbers. 
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Figure 2: Lamda Shock extracted from ONERA wing. 


downstream, along the top surface of the cylinder another shock forms due to the high angle 
of attack which causes the flow to reach supersonic speed on the leeward side of the body. 


Summary 

We have devised a method for shock extraction that uses knowledge of the flow to reduce 
computation time and noise. The basic steps of the method are as follows: 

1. Calculate the Mach component in the direction of the pressure gradient (Mj_). 

2. Verify that the region is undergoing compression. 

3. Verify that the pressure gradient within a region is consistent with the existence of a 
shock. 

4. Use the marching cube algorithm to search for Mi = 1.0. 

5. Disregard all surfaces whose normal is not in the same direction as the pressure gradi- 
ent. 


Future Work 

We plan on extending this work to other significant flow features. The extraction of vortices 
within a flow will be the first area of study followed by viewing skin friction along a body. 
The skin friction will be depicted through a simulation oil flow combined with pressure 
sensitive paint. 
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Figure 3: Bow shock with second shock on a hemi-spherical cylinder 
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Abstract. The increasing usage of computers in experimental and theoretical re- 
search has presented new challenges to the scientific visualization community. This 
is especially valid for the case of multi-dimensional datasets as vector fields which 
are produced in the field of chemical engineering when flow and concentration fields 
or heat transfer are calculated or measured. Vector fields can be visualized in detail 
using the “Line Integral Convolution” (LIC) Algorithm from Cabral and Leedom 
[1]. The focus of this paper is a discusion of basic principles and features of LIC, its 
usage on curvilinear grids and a comparison to traditional visualization techniques 
as vector arrows and stre aml ines. 


1 Introduction 

Traditional approaches for flow visualization are symbolic representations 
such as streamlines or vector arrows to visualize the direction of a flow field. 
In recent years several new techniques for vector field visualization have been 
introduced among the visualization with textures appears to be one of the 
most promising. There are several definitions for textures in literature. Com- 
monly Textures for flow visualization purposes are defined as surface struc- 
tures containing directional information. An example of a texture is shown in 
Fig. 1. For the first time van Wijk [2] used a texture to visualize the flow on a 
ship hull. With the introduction of the “Line Integral Convolution” -algorithm 
by Cabral and Leedom in 1993 [1] the generation of textures became much 
easier. Since 1993 there are published a number of papers concerning the 
application and extensions of the LIC-algorithm. Forssell and Cohen [3] ex- 
tended the LIC-algorithm to its application on curvilinear grids in 3D. As 
textures are by nature bidirectional, Delmarcelle and Hesselink [4,5] applied 
textures for the visualization of eigenvector fields of tensor fields. In order to 
visualize additionally the forward and backward direction and the velocity 
of a flow animation of textures has been introduced by Cabral and Leedom 
[1] and extended by several other research groups [3,6-9]. The generation of 
textures requires an input noise (cf. section 2). With a local variation of the 
frequency of the input noise Kiu and Banks [10] achieved special visual effects 
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Fig. 1. Flow field past a zylinder visualized with a texture 


using multi-frequency noise. Battke, Stalling and Hege [11] introduced an al- 
gorithm to calculate textures on arbitrary surfaces in 3D. In order to reduce 
the cost of computation for the calculation of textures, Stalling and Hege [6] 
introduced the algorithm "fastLIC 7- , which allows a more efficient calculation 
of textures. The visualization of unsteady flow fields with textures has been 
investigated by Forssell and Cohen [3] and by Shen and Kao [12]. 


2 Generation of Textures with Line Integral 
Convolution 

The basic principle of texture generation with the LIC-algorithm is to smear 
an input noise along the integral lines of a vector field. The input noise con- 
sists of random grey values for each pixel of the later generated texture. The 
smearing of the grey values can mathematically be expressed as a convolu- 
tion of the grey values along the integral lines of the vector field. If is the 
intensity- of the input noise (grey values) and c{$) ist the path or integral 
curve 

do-(s) v m 

as it?r w 

where the tangent of <7 has the direction of the vector field, then the intensity 
(grey) values I out of a texture can be calculated with the equation 

so — L 

I<mt(x,y) = J Iin{<r(s))h(s - s 0 )ds, with (x,y) = cr(s 0 )- (2) 

sq — L 
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In (2) h is the convolution filter function. This function has a strong influence 
on the appearance of a texture. The simplest filter function is a box filter with 


where L is the filter length in each direction from the output pixel. Thus 
/ = 2 * L is the whole filter length. The longer the filter length, the longer 
the black and white streaks will appear in the texture. If a box filter is used 
for the convolution (2) can be simplified to the summation 


where the intensity values of the pixels are only weighted by the path As 
the streamline passes through this pixel. In (4) (i.j) are the pixel coordinates 
at (x, y ) and P is a pixel of the set r of all pixels which are intersecting with 
the path a of the local streamline. 



Fig. 2. Visualization of a boundary layer flow over a fiat plate, (a) shows the multi- 
frequency input noise, (b) shows the generated texture 


The whole procedure to generate textures can be subdivided into four 


1. Calculation of the input noise, which consists of random intensity (grey) 
levels in the area, where the texture has to be generated. 

2. For each pixel a local streamline is calculated - starting at the center of 
this pixel in both directions. 

3. The intensity values of the input noise are integrated along the local 
streamline within the distance [— L. +L] and normalized by the integra- 
tion length. The pixel values will be weighted by a convolution filter 
function. 



( 3 ) 


Iout{x,y) — Iout{i:j) — ^ ] Iin(.P)ASi 

PCr(cr) 


( 4 ) 


steps: 
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4. The integration of the pixel values leads almost to an equalization of the 
intensity values. Therefore the contrast of the texture has to be increased. 

The result of this procedure is a black and white texture as it is depicted in 
Fig. 1. After that the texture can be colored in order to visualize an additional 
dimension of information (cf. Sect. 3). 

The type and frequeny of the input noise has a strong influence on the 
appearance of a texture. There are several algorithms to calculate the input 
noise [4,7]. The frequency of the noise can be locally varied in order to achieve 
an additional visual effect. In Fig. 2 the boundary layer flow on a flat plate 
is visualized using multi-frequency noise [10]. The velocity in the horizontal 
direction v x is a function of the vertical distance y to the plate as v x oc ^Jy. 
Besides the direction of the flow the texture in Fig. 2 gives an impression of 
the magitute of the velocity. 


Fig. 3. Visualization of the flow field of a viscoelatic fluid between two rotating 
rolls. The color encodes the pressure. Image by courtesy of Karsten Riest [13] 


3 Coloring Textures in Chemical Engineering 

With colored textures it is possible to visualize one more dimension of in- 
formation in the same image. The color can encode e.g. the magnitute of 
the velocity, the pressure, the volume fractions of the phases etc. In order to 
generate colored textures one first have to define a color map which gives for 


pressure: 

high 


low 
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volume fraction 
of the liquid: 

1.0 

0.5 

0.0 


Fig. 4. Flow field in a two-phase separator. The color encodes the volume fraction 
of the liquid. Image by courtesy of Matthias Creutz [14] 


each scalar value to be visualized a color, except black and white. Then the 
color values will be multiplied with the intensity value of the precalculated 
black and white texture [7]. In Fig. 3 the flow field between two rotating 
rolls is visualized where the color encodes the pressure. In such images the 
dependencies between pressure and flow topology can be studied easily. 

In the example in Fig. 4 the flow field in a two-phase separator is depicted. 
The flow is caculated with a commercial CFD programm. With the color 
map the volume fraction of the liquid phase is encoded. The texture shows 
the results on the plane in the middle of one section of the impeller. 

4 Textures on Curvilinear Surfaces 

Many grids for numerical calculations are curvilinear grids. It is possible to 
generate LIC-textures on curvilinear surfaces of curvilinear grids [3,7]. In 
order to generate a texture on a curvilinear surface, first the vector field on 
this surface is mathematically transformed into a 2D surface. We denote this 
surface the computational space. The transformation into the computational 
space can be performed using the Jacobian matrix. In the computational 
space the texture will be generated according to the algorithm described 
above. The resulting 2D-texture in computational space will be mapped back 
to the curvilinear surface in 3D. This procedure is illustrated on the left side 
in Fig. 5. In order to generate a texture of the flow on the first slice of the 
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multi-frequency noise 




LIC texture in 2D optimized LIC texture in 2D 

computational space computational space 


standard texture optimized texture 

mapped on surface mapped on surface 


Fig. 5. LIC-textures on curvilinear grids. Left: Mapping a 2D texture on a curvi- 
linear surface in 3D. Right: Optimize the streak size of the 3D-texture using multi- 
frequency noise 
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numerical grid of a hemishere cylinder (top), a LIC-texture in computational 
space is calculated (middle) and then mapped back to the curvilinear surface 
(bottom). 

This procedure works well if the grid density is similar in all areas of 
the curvilinear surface. But for the hemisphere cylinder the grid density is 
much higher at the top of the hemishere than at the end of the zylinder. 
Thus the streaks which are homogeneous in the computational space vary 
significantly on the curvilinear surface in 3D (Fig. 5 bottom left). This can 
be compensated by using multi-frequency noise. The frequency of the noise 
is calculated high, where the grid density is low and vice versa (Fig. 5 top 
right). Thus locally different streak widths are calculated for the texture in 
computational space (Fig. 5 middle left). If this texture is mapped onto the 
curvilinear surface it results in an almost uniform streak width in 3D (Fig. 5 
bottom right). 


Table 1 . Advantages and disadvantages of flow visualization with textures com- 
pared to traditional visualization techniques 


visualization 

technique 

advantages 

disadvantages 

vector arrows 

simple depiction 
direction clear 

length of arrow shows mag- 
nitude of velocity 

details difficult to visualize 
discretisation necessary 

visual clutter for high den- 
sity of arrows 

streamlines 

visualization of flow path 

quantitative comparisons of 
different flow fields possible 

details difficult to visualize 

discretisation necessary, 

problems for convergence / 
divergence 

bidirectional 

textures 

visual appealing 

highest possible resolution of 
flow fields, details depicted 

good representation of flow 
topology- 

continuous representation 

high computational cost 
bidirectional 


5 Textures versus Traditional Visualization Techniques 

The advantages and disadvantages of textures compared to traditional flow 
visualization techniques, among which vector arrows and streamlines are the 
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most common, are shown in Table 4. One of the major strength of textures is 
that with pixel-oriented textures the highest possible resolution of a flow field 
can be achieved. Thus even details like secondary vortices can be resolved and 
the flow topology is easy to extract graphically. The disadvantage of textures 
for flow visualization is their high computational cost and the bidirectional 
nature of textures. The later disadvantage can be overcome by animation 
of textures [6,7]. The bidirectional nature of textures makes them ideal for 
the visualization of bidirectional vector fields like eigenvector fields of tensor 
fields [5]. 
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